Fixes for coverity checks.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliCentralMCCorrectionsTask.cxx
CommitLineData
600b313f 1//
2// Calculate the corrections in the central regions
3//
4// Inputs:
5// - AliESDEvent
6//
7// Outputs:
8// - AliAODCentralMult
9//
10// Histograms
11//
12// Corrections used
13//
14#include "AliCentralMCCorrectionsTask.h"
15#include "AliTriggerAnalysis.h"
16#include "AliPhysicsSelection.h"
17#include "AliLog.h"
18#include "AliHeader.h"
19#include "AliGenEventHeader.h"
20#include "AliESDEvent.h"
21#include "AliAODHandler.h"
22#include "AliMultiplicity.h"
23#include "AliInputEventHandler.h"
24#include "AliStack.h"
25#include "AliMCEvent.h"
26#include "AliAODForwardMult.h"
27#include "AliCentralCorrSecondaryMap.h"
28#include <TH1.h>
29#include <TH2D.h>
30#include <TDirectory.h>
31#include <TList.h>
32#include <TROOT.h>
33#include <iostream>
34
35//====================================================================
36namespace {
37 const char* GetEventName(Bool_t tr, Bool_t vtx)
38 {
39 return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : ""));
40 }
41 const char* GetHitsName(UShort_t d, Char_t r)
42 {
43 return Form("hitsSPD%d%c", d, r);
44 }
45 const char* GetStripsName(UShort_t d, Char_t r)
46 {
47 return Form("stripsSPD%d%c", d, r);
48 }
49 const char* GetPrimaryName(Char_t r, Bool_t trVtx)
50 {
51 return Form("primaries%s%s",
52 ((r == 'I' || r == 'i') ? "Inner" : "Outer"),
53 (trVtx ? "TrVtx" : "All"));
54 }
55}
56
57//====================================================================
58AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask()
59 : AliAnalysisTaskSE(),
60 fInspector(),
61 fTrackDensity(),
62 fVtxBins(0),
63 fFirstEvent(true),
64 fHEvents(0),
65 fHEventsTr(0),
66 fHEventsTrVtx(0),
67 fVtxAxis(),
68 fEtaAxis(),
69 fList(),
4a9f6ae6 70 fNPhiBins(20)
600b313f 71{
72 //
73 // Constructor
74 //
75 // Parameters:
76 // name Name of task
77 //
78}
79
80//____________________________________________________________________
81AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const char* name)
82 : AliAnalysisTaskSE(name),
83 fInspector("eventInspector"),
84 fTrackDensity("trackDensity"),
85 fVtxBins(0),
86 fFirstEvent(true),
87 fHEvents(0),
88 fHEventsTr(0),
89 fHEventsTrVtx(0),
90 fVtxAxis(10,-10,10),
91 fEtaAxis(200,-4,6),
92 fList(),
4a9f6ae6 93 fNPhiBins(20)
600b313f 94{
95 //
96 // Constructor
97 //
98 // Parameters:
99 // name Name of task
100 //
101 DefineOutput(1, TList::Class());
102 DefineOutput(2, TList::Class());
103}
104
105//____________________________________________________________________
106AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const AliCentralMCCorrectionsTask& o)
107 : AliAnalysisTaskSE(o),
108 fInspector(o.fInspector),
109 fTrackDensity(),
110 fVtxBins(0),
111 fFirstEvent(o.fFirstEvent),
112 fHEvents(o.fHEvents),
113 fHEventsTr(o.fHEventsTr),
114 fHEventsTrVtx(o.fHEventsTrVtx),
115 fVtxAxis(10,-10,10),
116 fEtaAxis(200,-4,6),
117 fList(o.fList),
118 fNPhiBins(o.fNPhiBins)
119{
120 //
121 // Copy constructor
122 //
123 // Parameters:
124 // o Object to copy from
125 //
126}
127
128//____________________________________________________________________
129AliCentralMCCorrectionsTask&
130AliCentralMCCorrectionsTask::operator=(const AliCentralMCCorrectionsTask& o)
131{
132 //
133 // Assignment operator
134 //
135 // Parameters:
136 // o Object to assign from
137 //
138 // Return:
139 // Reference to this object
140 //
d015ecfe 141 if (&o == this) return *this;
600b313f 142 fInspector = o.fInspector;
143 fTrackDensity = o.fTrackDensity;
144 fVtxBins = o.fVtxBins;
145 fFirstEvent = o.fFirstEvent;
146 fHEvents = o.fHEvents;
147 fHEventsTr = o.fHEventsTr;
148 fHEventsTrVtx = o.fHEventsTrVtx;
149 SetVertexAxis(o.fVtxAxis);
150 SetEtaAxis(o.fEtaAxis);
151 fNPhiBins = o.fNPhiBins;
152
153 return *this;
154}
155
156//____________________________________________________________________
157void
158AliCentralMCCorrectionsTask::Init()
159{
160 //
161 // Initialize the task
162 //
163 //
164}
165
166//____________________________________________________________________
167void
168AliCentralMCCorrectionsTask::SetVertexAxis(Int_t nBin, Double_t min,
169 Double_t max)
170{
171 //
172 // Set the vertex axis to use
173 //
174 // Parameters:
175 // nBins Number of bins
176 // vzMin Least @f$z@f$ coordinate of interation point
177 // vzMax Largest @f$z@f$ coordinate of interation point
178 //
179 if (max < min) max = -min;
180 if (min < max) {
181 Double_t tmp = min;
182 min = max;
183 max = tmp;
184 }
185 if (min < -10)
186 AliWarning(Form("Minimum vertex %f < -10, make sure you want this",min));
187 if (max > +10)
188 AliWarning(Form("Minimum vertex %f > +10, make sure you want this",max));
189 fVtxAxis.Set(nBin, min, max);
190}
191//____________________________________________________________________
192void
193AliCentralMCCorrectionsTask::SetVertexAxis(const TAxis& axis)
194{
195 //
196 // Set the vertex axis to use
197 //
198 // Parameters:
199 // axis Axis
200 //
201 SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
202}
203
204//____________________________________________________________________
205void
206AliCentralMCCorrectionsTask::SetEtaAxis(Int_t nBin, Double_t min, Double_t max)
207{
208 //
209 // Set the eta axis to use
210 //
211 // Parameters:
212 // nBins Number of bins
213 // vzMin Least @f$\eta@f$
214 // vzMax Largest @f$\eta@f$
215 //
216 if (max < min) max = -min;
217 if (min < max) {
218 Double_t tmp = min;
219 min = max;
220 max = tmp;
221 }
222 if (min < -4)
223 AliWarning(Form("Minimum eta %f < -4, make sure you want this",min));
224 if (max > +6)
225 AliWarning(Form("Minimum eta %f > +6, make sure you want this",max));
226 fVtxAxis.Set(nBin, min, max);
227}
228//____________________________________________________________________
229void
230AliCentralMCCorrectionsTask::SetEtaAxis(const TAxis& axis)
231{
232 //
233 // Set the eta axis to use
234 //
235 // Parameters:
236 // axis Axis
237 //
238 SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
239}
240
241//____________________________________________________________________
242void
243AliCentralMCCorrectionsTask::DefineBins(TList* l)
244{
245 if (!fVtxBins) fVtxBins = new TObjArray(fVtxAxis.GetNbins(), 1);
246 if (fVtxBins->GetEntries() > 0) return;
247
248 fVtxBins->SetOwner();
249 for (Int_t i = 1; i <= fVtxAxis.GetNbins(); i++) {
250 Double_t low = fVtxAxis.GetBinLowEdge(i);
251 Double_t high = fVtxAxis.GetBinUpEdge(i);
252 VtxBin* bin = new VtxBin(low, high, fEtaAxis, fNPhiBins);
253 fVtxBins->AddAt(bin, i);
254 bin->DefineOutput(l);
255 }
256}
257
258//____________________________________________________________________
259void
260AliCentralMCCorrectionsTask::UserCreateOutputObjects()
261{
262 //
263 // Create output objects
264 //
265 //
266 fList = new TList;
267 fList->SetOwner();
268 fList->SetName(Form("%sSums", GetName()));
269
270 DefineBins(fList);
271
272 fHEvents = new TH1I(GetEventName(false,false),
273 "Number of all events",
274 fVtxAxis.GetNbins(),
275 fVtxAxis.GetXmin(),
276 fVtxAxis.GetXmax());
277 fHEvents->SetXTitle("v_{z} [cm]");
278 fHEvents->SetYTitle("# of events");
279 fHEvents->SetFillColor(kBlue+1);
280 fHEvents->SetFillStyle(3001);
281 fHEvents->SetDirectory(0);
282 fList->Add(fHEvents);
283
284 fHEventsTr = new TH1I(GetEventName(true, false),
285 "Number of triggered events",
286 fVtxAxis.GetNbins(),
287 fVtxAxis.GetXmin(),
288 fVtxAxis.GetXmax());
289 fHEventsTr->SetXTitle("v_{z} [cm]");
290 fHEventsTr->SetYTitle("# of events");
291 fHEventsTr->SetFillColor(kRed+1);
292 fHEventsTr->SetFillStyle(3001);
293 fHEventsTr->SetDirectory(0);
294 fList->Add(fHEventsTr);
295
296 fHEventsTrVtx = new TH1I(GetEventName(true,true),
297 "Number of events w/trigger and vertex",
298 fVtxAxis.GetNbins(),
299 fVtxAxis.GetXmin(),
300 fVtxAxis.GetXmax());
301 fHEventsTrVtx->SetXTitle("v_{z} [cm]");
302 fHEventsTrVtx->SetYTitle("# of events");
303 fHEventsTrVtx->SetFillColor(kBlue+1);
304 fHEventsTrVtx->SetFillStyle(3001);
305 fHEventsTrVtx->SetDirectory(0);
306 fList->Add(fHEventsTrVtx);
307
308 // Copy axis objects to output
309 TH1* vtxAxis = new TH1D("vtxAxis", "Vertex axis",
310 fVtxAxis.GetNbins(),
311 fVtxAxis.GetXmin(),
312 fVtxAxis.GetXmax());
313 TH1* etaAxis = new TH1D("etaAxis", "Eta axis",
314 fEtaAxis.GetNbins(),
315 fEtaAxis.GetXmin(),
316 fEtaAxis.GetXmax());
317 fList->Add(vtxAxis);
318 fList->Add(etaAxis);
319
320 AliInfo(Form("Initialising sub-routines: %p, %p",
321 &fInspector, &fTrackDensity));
322 fInspector.DefineOutput(fList);
323 fInspector.Init(fVtxAxis);
324 fTrackDensity.DefineOutput(fList);
325
326 PostData(1, fList);
327}
328//____________________________________________________________________
329void
330AliCentralMCCorrectionsTask::UserExec(Option_t*)
331{
332 //
333 // Process each event
334 //
335 // Parameters:
336 // option Not used
337 //
338
339 // Get the input data - MC event
340 AliMCEvent* mcEvent = MCEvent();
341 if (!mcEvent) {
342 AliWarning("No MC event found");
343 return;
344 }
345
346 // Get the input data - ESD event
347 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
348 if (!esd) {
349 AliWarning("No ESD event found for input event");
350 return;
351 }
352
353 //--- Read run information -----------------------------------------
354 if (fFirstEvent && esd->GetESDRun()) {
355 fInspector.ReadRunDetails(esd);
356
357 AliInfo(Form("Initializing with parameters from the ESD:\n"
358 " AliESDEvent::GetBeamEnergy() ->%f\n"
359 " AliESDEvent::GetBeamType() ->%s\n"
360 " AliESDEvent::GetCurrentL3() ->%f\n"
361 " AliESDEvent::GetMagneticField()->%f\n"
362 " AliESDEvent::GetRunNumber() ->%d\n",
363 esd->GetBeamEnergy(),
364 esd->GetBeamType(),
365 esd->GetCurrentL3(),
366 esd->GetMagneticField(),
367 esd->GetRunNumber()));
368
369 Print();
370 fFirstEvent = false;
371 }
372
373 // Some variables
374 UInt_t triggers; // Trigger bits
375 Bool_t lowFlux; // Low flux flag
376 UShort_t iVz; // Vertex bin from ESD
377 Double_t vZ; // Z coordinate from ESD
378 Double_t cent; // Centrality
379 UShort_t iVzMc; // Vertex bin from MC
380 Double_t vZMc; // Z coordinate of IP vertex from MC
381 Double_t b; // Impact parameter
382 Int_t nPart; // Number of participants
383 Int_t nBin; // Number of binary collisions
384 Double_t phiR; // Reaction plane from MC
385 UShort_t nClusters;// Number of SPD clusters
386 // Process the data
387 UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ,
388 cent, nClusters);
389 fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc,
390 b, nPart, nBin, phiR);
391
392 Bool_t isInel = triggers & AliAODForwardMult::kInel;
393 Bool_t hasVtx = retESD == AliFMDMCEventInspector::kOk;
394
395 // Fill the event count histograms
396 if (isInel) fHEventsTr->Fill(vZMc);
397 if (isInel && hasVtx) fHEventsTrVtx->Fill(vZMc);
398 fHEvents->Fill(vZMc);
399
400 // Now find our vertex bin object
401 VtxBin* bin = 0;
402 if (iVzMc > 0 && iVzMc <= fVtxAxis.GetNbins())
403 bin = static_cast<VtxBin*>(fVtxBins->At(iVzMc));
404 if (!bin) {
405 // AliError(Form("No vertex bin object @ %d (%f)", iVzMc, vZMc));
406 return;
407 }
408
409 // Now process our input data and store in own ESD object
410 fTrackDensity.Calculate(*mcEvent, vZMc, *bin->fHits, bin->fPrimary);
411 bin->fCounts->Fill(0.5);
412}
413
414//____________________________________________________________________
415void
416AliCentralMCCorrectionsTask::Terminate(Option_t*)
417{
418 //
419 // End of job
420 //
421 // Parameters:
422 // option Not used
423 //
424
425 fList = dynamic_cast<TList*>(GetOutputData(1));
426 if (!fList) {
427 AliError("No output list defined");
428 return;
429 }
430
431 // Output list
432 TList* output = new TList;
433 output->SetOwner();
434 output->SetName(Form("%sResults", GetName()));
435
436 // --- Fill correction object --------------------------------------
437 AliCentralCorrSecondaryMap* corr = new AliCentralCorrSecondaryMap;
438 corr->SetVertexAxis(fVtxAxis);
439 // corr->SetEtaAxis(fEtaAxis);
440
441 TIter next(fVtxBins);
442 VtxBin* bin = 0;
443 UShort_t iVz = 1;
444 while ((bin = static_cast<VtxBin*>(next())))
445 bin->Finish(fList, output, iVz++, corr);
446
447 output->Add(corr);
448
449 PostData(2, output);
450}
451
452//____________________________________________________________________
453void
454AliCentralMCCorrectionsTask::Print(Option_t* option) const
455{
456 std::cout << ClassName() << "\n"
457 << " Vertex bins: " << fVtxAxis.GetNbins() << '\n'
458 << " Vertex range: [" << fVtxAxis.GetXmin()
459 << "," << fVtxAxis.GetXmax() << "]\n"
460 << " Eta bins: " << fEtaAxis.GetNbins() << '\n'
461 << " Eta range: [" << fEtaAxis.GetXmin()
462 << "," << fEtaAxis.GetXmax() << "]\n"
463 << " # of phi bins: " << fNPhiBins
464 << std::endl;
465 gROOT->IncreaseDirLevel();
466 fInspector.Print(option);
467 fTrackDensity.Print(option);
468 gROOT->DecreaseDirLevel();
469}
470
471//====================================================================
472const char*
473AliCentralMCCorrectionsTask::VtxBin::BinName(Double_t low,
474 Double_t high)
475{
476 static TString buf;
477 buf = Form("vtx%+05.1f_%+05.1f", low, high);
478 buf.ReplaceAll("+", "p");
479 buf.ReplaceAll("-", "m");
480 buf.ReplaceAll(".", "d");
481 return buf.Data();
482}
483
484
485//____________________________________________________________________
486AliCentralMCCorrectionsTask::VtxBin::VtxBin()
487 : fHits(0),
488 fPrimary(0),
489 fCounts(0)
490{
491}
492//____________________________________________________________________
493AliCentralMCCorrectionsTask::VtxBin::VtxBin(Double_t low,
494 Double_t high,
495 const TAxis& axis,
496 UShort_t nPhi)
497 : TNamed(BinName(low, high),
498 Form("%+5.1fcm<v_{z}<%+5.1fcm", low, high)),
499 fHits(0),
500 fPrimary(0),
501 fCounts(0)
502{
503 fPrimary = new TH2D("primary", "Primaries",
504 axis.GetNbins(), axis.GetXmin(), axis.GetXmax(),
505 nPhi, 0, 2*TMath::Pi());
506 fPrimary->SetXTitle("#eta");
507 fPrimary->SetYTitle("#varphi [radians]");
508 fPrimary->Sumw2();
509 fPrimary->SetDirectory(0);
510
511 fHits = static_cast<TH2D*>(fPrimary->Clone("hits"));
512 fHits->SetTitle("Hits");
513 fHits->SetDirectory(0);
514
515 fCounts = new TH1D("counts", "Counts", 1, 0, 1);
516 fCounts->SetXTitle("Events");
517 fCounts->SetYTitle("# of Events");
518 fCounts->SetDirectory(0);
519}
520
521//____________________________________________________________________
522AliCentralMCCorrectionsTask::VtxBin::VtxBin(const VtxBin& o)
523 : TNamed(o),
524 fHits(0),
525 fPrimary(0),
526 fCounts(0)
527{
528 if (o.fHits) {
529 fHits = static_cast<TH2D*>(o.fHits->Clone());
530 fHits->SetDirectory(0);
531 }
532 if (o.fPrimary) {
533 fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
534 fPrimary->SetDirectory(0);
535 }
536 if (o.fCounts) {
537 fCounts = static_cast<TH1D*>(o.fCounts->Clone());
538 fCounts->SetDirectory(0);
539 }
540}
541
542//____________________________________________________________________
543AliCentralMCCorrectionsTask::VtxBin&
544AliCentralMCCorrectionsTask::VtxBin::operator=(const VtxBin& o)
545{
d015ecfe 546 if (&o == this) return *this;
600b313f 547 TNamed::operator=(o);
548 fHits = 0;
549 fPrimary = 0;
550 fCounts = 0;
551 if (o.fHits) {
552 fHits = static_cast<TH2D*>(o.fHits->Clone());
553 fHits->SetDirectory(0);
554 }
555 if (o.fPrimary) {
556 fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
557 fPrimary->SetDirectory(0);
558 }
559 if (o.fCounts) {
560 fCounts = static_cast<TH1D*>(o.fCounts->Clone());
561 fCounts->SetDirectory(0);
562 }
563 return *this;
564}
565
566//____________________________________________________________________
567void
568AliCentralMCCorrectionsTask::VtxBin::DefineOutput(TList* l)
569{
570 TList* d = new TList;
571 d->SetName(GetName());
572 d->SetOwner();
573 l->Add(d);
574
575 d->Add(fHits);
576 d->Add(fPrimary);
577 d->Add(fCounts);
578}
579//____________________________________________________________________
580void
581AliCentralMCCorrectionsTask::VtxBin::Finish(const TList* input,
582 TList* output,
583 UShort_t iVz,
584 AliCentralCorrSecondaryMap* map)
585{
586 TList* out = new TList;
587 out->SetName(GetName());
588 out->SetOwner();
589 output->Add(out);
590
591 TList* l = static_cast<TList*>(input->FindObject(GetName()));
592 if (!l) {
593 AliError(Form("List %s not found in %s", GetName(), input->GetName()));
594 return;
595 }
596
597
598 TH2D* hits = static_cast<TH2D*>(l->FindObject("hits"));
599 TH2D* prim = static_cast<TH2D*>(l->FindObject("primary"));
600 if (!hits || !prim) {
601 AliError(Form("Missing histograms: %p, %p", hits, prim));
602 return;
603 }
604
605 TH2D* h = static_cast<TH2D*>(hits->Clone("bgCorr"));
606 h->SetDirectory(0);
984a7263 607 h->Divide(prim);
600b313f 608
609 map->SetCorrection(iVz, h);
610
611 out->Add(h);
612}
613
614//
615// EOF
616//