More code clean up.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardMultiplicityTask.cxx
CommitLineData
0bd4b00f 1#include "AliForwardMultiplicityTask.h"
7e4038b5 2#include "AliTriggerAnalysis.h"
3#include "AliPhysicsSelection.h"
4#include "AliLog.h"
0bd4b00f 5// #include "AliFMDAnaParameters.h"
7e4038b5 6#include "AliESDEvent.h"
7#include "AliAODHandler.h"
8#include "AliMultiplicity.h"
9#include "AliInputEventHandler.h"
0bd4b00f 10#include "AliForwardCorrectionManager.h"
11#include "AliAnalysisManager.h"
7e4038b5 12#include <TH1.h>
13#include <TDirectory.h>
14#include <TTree.h>
0bd4b00f 15#include <TROOT.h>
16#include <iostream>
17#include <iomanip>
7e4038b5 18
19//====================================================================
0bd4b00f 20AliForwardMultiplicityTask::AliForwardMultiplicityTask()
7e4038b5 21 : AliAnalysisTaskSE(),
0bd4b00f 22 fEnableLowFlux(true),
7e4038b5 23 fHData(0),
24 fFirstEvent(true),
7e4038b5 25 fESDFMD(),
26 fHistos(),
27 fAODFMD(),
fffea31d 28 fEventInspector(),
29 fEnergyFitter(),
7e4038b5 30 fSharingFilter(),
31 fDensityCalculator(),
32 fCorrections(),
33 fHistCollector(),
fffea31d 34 fList(0)
7e4038b5 35{
36}
37
38//____________________________________________________________________
0bd4b00f 39AliForwardMultiplicityTask::AliForwardMultiplicityTask(const char* name)
7e4038b5 40 : AliAnalysisTaskSE(name),
0bd4b00f 41 fEnableLowFlux(true),
7e4038b5 42 fHData(0),
43 fFirstEvent(true),
7e4038b5 44 fESDFMD(),
45 fHistos(),
46 fAODFMD(kTRUE),
fffea31d 47 fEventInspector("event"),
48 fEnergyFitter("energy"),
7e4038b5 49 fSharingFilter("sharing"),
50 fDensityCalculator("density"),
51 fCorrections("corrections"),
52 fHistCollector("collector"),
fffea31d 53 fList(0)
7e4038b5 54{
55 DefineOutput(1, TList::Class());
7e4038b5 56}
57
58//____________________________________________________________________
0bd4b00f 59AliForwardMultiplicityTask::AliForwardMultiplicityTask(const AliForwardMultiplicityTask& o)
7e4038b5 60 : AliAnalysisTaskSE(o),
0bd4b00f 61 fEnableLowFlux(o.fEnableLowFlux),
7e4038b5 62 fHData(o.fHData),
63 fFirstEvent(true),
7e4038b5 64 fESDFMD(o.fESDFMD),
65 fHistos(o.fHistos),
66 fAODFMD(o.fAODFMD),
fffea31d 67 fEventInspector(o.fEventInspector),
68 fEnergyFitter(o.fEnergyFitter),
7e4038b5 69 fSharingFilter(o.fSharingFilter),
70 fDensityCalculator(o.fDensityCalculator),
71 fCorrections(o.fCorrections),
72 fHistCollector(o.fHistCollector),
fffea31d 73 fList(o.fList)
7e4038b5 74{
0bd4b00f 75 DefineOutput(1, TList::Class());
7e4038b5 76}
77
78//____________________________________________________________________
0bd4b00f 79AliForwardMultiplicityTask&
80AliForwardMultiplicityTask::operator=(const AliForwardMultiplicityTask& o)
7e4038b5 81{
fffea31d 82 AliAnalysisTaskSE::operator=(o);
83
0bd4b00f 84 fEnableLowFlux = o.fEnableLowFlux;
7e4038b5 85 fHData = o.fHData;
86 fFirstEvent = o.fFirstEvent;
fffea31d 87 fEventInspector = o.fEventInspector;
88 fEnergyFitter = o.fEnergyFitter;
7e4038b5 89 fSharingFilter = o.fSharingFilter;
90 fDensityCalculator = o.fDensityCalculator;
91 fCorrections = o.fCorrections;
92 fHistCollector = o.fHistCollector;
93 fHistos = o.fHistos;
94 fAODFMD = o.fAODFMD;
95 fList = o.fList;
7e4038b5 96
97 return *this;
98}
99
fffea31d 100//____________________________________________________________________
101void
0bd4b00f 102AliForwardMultiplicityTask::SetDebug(Int_t dbg)
fffea31d 103{
104 fEventInspector.SetDebug(dbg);
105 fEnergyFitter.SetDebug(dbg);
106 fSharingFilter.SetDebug(dbg);
107 fDensityCalculator.SetDebug(dbg);
108 fCorrections.SetDebug(dbg);
109 fHistCollector.SetDebug(dbg);
110}
111
7e4038b5 112//____________________________________________________________________
113void
0bd4b00f 114AliForwardMultiplicityTask::Init()
7e4038b5 115{
116 fFirstEvent = true;
117}
118
119//____________________________________________________________________
120void
0bd4b00f 121AliForwardMultiplicityTask::InitializeSubs()
7e4038b5 122{
7e4038b5 123
0bd4b00f 124 // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
125 // pars->Init(kTRUE);
126
127 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
128 fcm.Init(fEventInspector.GetCollisionSystem(),
129 fEventInspector.GetEnergy(),
130 fEventInspector.GetField());
131 // Check that we have the energy loss fits, needed by
132 // AliFMDSharingFilter
133 // AliFMDDensityCalculator
134 if (!fcm.GetELossFit()) {
135 AliFatal(Form("No energy loss fits"));
136 return;
137 }
138 // Check that we have the double hit correction - (optionally) used by
139 // AliFMDDensityCalculator
140 if (!fcm.GetDoubleHit()) {
141 AliWarning("No double hit corrections");
142 }
143 // Check that we have the secondary maps, needed by
144 // AliFMDCorrections
145 // AliFMDHistCollector
146 if (!fcm.GetSecondaryMap()) {
147 AliFatal("No secondary corrections");
148 return;
149 }
150 // Check that we have the vertex bias correction, needed by
151 // AliFMDCorrections
152 if (!fcm.GetVertexBias()) {
153 AliFatal("No event vertex bias corrections");
154 return;
155 }
156 // Check that we have the merging efficiencies, optionally used by
157 // AliFMDCorrections
158 if (!fcm.GetMergingEfficiency()) {
159 AliWarning("No merging efficiencies");
160 }
161
162
163 const TAxis* pe = fcm.GetEtaAxis();
164 const TAxis* pv = fcm.GetVertexAxis();
165 if (!pe) AliFatal("No eta axis defined");
166 if (!pv) AliFatal("No vertex axis defined");
fffea31d 167
0bd4b00f 168 // TAxis e(pars->GetNetaBins(), pars->GetEtaMin(), pars->GetEtaMax());
169 // TAxis v(pars->GetNvtxBins(), -pars->GetVtxCutZ(), pars->GetVtxCutZ());
170 // pv->Dump();
171
172 fHistos.Init(*pe);
173 fAODFMD.Init(*pe);
7e4038b5 174
175 fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
176 fHData->SetStats(0);
177 fHData->SetDirectory(0);
9d99b0dd 178 fList->Add(fHData);
179
0bd4b00f 180 fEnergyFitter.Init(*pe);
181 fEventInspector.Init(*pv);
182 fHistCollector.Init(*pv);
183
184 this->Print();
7e4038b5 185}
186
187//____________________________________________________________________
188void
0bd4b00f 189AliForwardMultiplicityTask::UserCreateOutputObjects()
7e4038b5 190{
191 fList = new TList;
192
193 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
0bd4b00f 194 AliAODHandler* ah =
195 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
fffea31d 196 if (!ah) AliFatal("No AOD output handler set in analysis manager");
7e4038b5 197
198
199 TObject* obj = &fAODFMD;
200 ah->AddBranch("AliAODForwardMult", &obj);
201
fffea31d 202 fEventInspector.DefineOutput(fList);
203 fEnergyFitter.DefineOutput(fList);
9d99b0dd 204 fSharingFilter.DefineOutput(fList);
205 fDensityCalculator.DefineOutput(fList);
206 fCorrections.DefineOutput(fList);
0bd4b00f 207
208 PostData(1, fList);
7e4038b5 209}
210//____________________________________________________________________
211void
0bd4b00f 212AliForwardMultiplicityTask::UserExec(Option_t*)
7e4038b5 213{
0bd4b00f 214 // static Int_t cnt = 0;
215 // cnt++;
7e4038b5 216 // Get the input data
217 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
0bd4b00f 218 // AliInfo(Form("Event # %6d (esd=%p)", cnt, esd));
7e4038b5 219 if (!esd) {
220 AliWarning("No ESD event found for input event");
221 return;
222 }
223
7e4038b5 224 // On the first event, initialize the parameters
fffea31d 225 if (fFirstEvent && esd->GetESDRun()) {
0bd4b00f 226 fEventInspector.ReadRunDetails(esd);
227
fffea31d 228 AliInfo(Form("Initializing with parameters from the ESD:\n"
229 " AliESDEvent::GetBeamEnergy() ->%f\n"
230 " AliESDEvent::GetBeamType() ->%s\n"
231 " AliESDEvent::GetCurrentL3() ->%f\n"
232 " AliESDEvent::GetMagneticField()->%f\n"
233 " AliESDEvent::GetRunNumber() ->%d\n",
234 esd->GetBeamEnergy(),
235 esd->GetBeamType(),
236 esd->GetCurrentL3(),
237 esd->GetMagneticField(),
238 esd->GetRunNumber()));
0bd4b00f 239
240
241
242 // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
243 // pars->SetParametersFromESD(esd);
244 // pars->PrintStatus();
7e4038b5 245 fFirstEvent = false;
246
247 InitializeSubs();
248 }
249 // Clear stuff
250 fHistos.Clear();
251 fESDFMD.Clear();
252 fAODFMD.Clear();
253
fffea31d 254 Bool_t lowFlux = kFALSE;
255 UInt_t triggers = 0;
0bd4b00f 256 UShort_t ivz = 0;
fffea31d 257 Double_t vz = 0;
258 UInt_t found = fEventInspector.Process(esd, triggers, lowFlux, ivz, vz);
259 if (found & AliFMDEventInspector::kNoEvent) return;
260 if (found & AliFMDEventInspector::kNoTriggers) return;
261
262 // Set trigger bits, and mark this event for storage
9d99b0dd 263 fAODFMD.SetTriggerBits(triggers);
7e4038b5 264 MarkEventForStore();
265
fffea31d 266 if (found & AliFMDEventInspector::kNoSPD) return;
267 if (found & AliFMDEventInspector::kNoFMD) return;
268 if (found & AliFMDEventInspector::kNoVertex) return;
269 fAODFMD.SetIpZ(vz);
7e4038b5 270
fffea31d 271 if (found & AliFMDEventInspector::kBadVertex) return;
7e4038b5 272
0bd4b00f 273 // We we do not want to use low flux specific code, we disable it here.
274 if (!fEnableLowFlux) lowFlux = false;
275
fffea31d 276 // Get FMD data
277 AliESDFMD* esdFMD = esd->GetFMDData();
278 // Apply the sharing filter (or hit merging or clustering if you like)
279 if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) {
280 AliWarning("Sharing filter failed!");
7e4038b5 281 return;
282 }
283
fffea31d 284 // Do the energy stuff
285 if (!fEnergyFitter.Accumulate(*esdFMD, triggers & AliAODForwardMult::kEmpty)){
286 AliWarning("Energy fitter failed");
7e4038b5 287 return;
288 }
289
290 // Calculate the inclusive charged particle density
291 if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
292 AliWarning("Density calculator failed!");
293 return;
294 }
295
296 // Do the secondary and other corrections.
297 if (!fCorrections.Correct(fHistos, ivz)) {
298 AliWarning("Corrections failed");
299 return;
300 }
301
302 if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
303 AliWarning("Histogram collector failed");
304 return;
305 }
306
307 if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
308 fHData->Add(&(fAODFMD.GetHistogram()));
2d68d438 309
310 PostData(1, fList);
7e4038b5 311}
312
313//____________________________________________________________________
314void
0bd4b00f 315AliForwardMultiplicityTask::Terminate(Option_t*)
7e4038b5 316{
317 TList* list = dynamic_cast<TList*>(GetOutputData(1));
318 if (!list) {
2d68d438 319 AliError(Form("No output list defined (%p)", GetOutputData(1)));
320 if (GetOutputData(1)) GetOutputData(1)->Print();
7e4038b5 321 return;
322 }
9d99b0dd 323
324 // Get our histograms from the container
fffea31d 325 TH1I* hEventsTr = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
326 TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
327 TH1I* hTriggers = 0;
328 if (!fEventInspector.FetchHistograms(list, hEventsTr,
329 hEventsTrVtx, hTriggers)) {
330 AliError(Form("Didn't get histograms from event selector "
331 "(hEventsTr=%p,hEventsTrVtx=%p)",
332 hEventsTr, hEventsTrVtx));
333 list->ls();
334 return;
335 }
336
9d99b0dd 337 TH2D* hData = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
fffea31d 338 if (!hData) {
339 AliError(Form("Couldn't get our summed histogram from output "
340 "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
2d68d438 341 list->ls();
342 return;
343 }
9d99b0dd 344
7e4038b5 345 // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
9d99b0dd 346 TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
0bd4b00f 347 TH1D* norm = hData->ProjectionX("norm", 0, 1, "");
7e4038b5 348 dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
349 dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
350 dNdeta->Divide(norm);
351 dNdeta->SetStats(0);
9d99b0dd 352 dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
7e4038b5 353 "width");
7e4038b5 354 list->Add(dNdeta);
0bd4b00f 355 list->Add(norm);
fffea31d 356
357 fEnergyFitter.Fit(list);
9d99b0dd 358 fSharingFilter.ScaleHistograms(list,hEventsTr->Integral());
359 fDensityCalculator.ScaleHistograms(list,hEventsTrVtx->Integral());
360 fCorrections.ScaleHistograms(list,hEventsTrVtx->Integral());
7e4038b5 361}
362
363//____________________________________________________________________
364void
0bd4b00f 365AliForwardMultiplicityTask::MarkEventForStore() const
7e4038b5 366{
367 // Make sure the AOD tree is filled
368 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
369 AliAODHandler* ah =
370 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
371 if (!ah)
372 AliFatal("No AOD output handler set in analysis manager");
373
374 ah->SetFillAOD(kTRUE);
375}
7e4038b5 376
377//____________________________________________________________________
378void
0bd4b00f 379AliForwardMultiplicityTask::Print(Option_t* option) const
7e4038b5 380{
0bd4b00f 381 std::cout << "AliForwardMultiplicityTask: " << GetName() << "\n"
382 << " Enable low flux code: " << (fEnableLowFlux ? "yes" : "no")
383 << std::endl;
384 gROOT->IncreaseDirLevel();
385 fEventInspector .Print(option);
386 fEnergyFitter .Print(option);
387 fSharingFilter .Print(option);
388 fDensityCalculator.Print(option);
389 fCorrections .Print(option);
390 fHistCollector .Print(option);
391 gROOT->DecreaseDirLevel();
7e4038b5 392}
393
394//
395// EOF
396//