New MC correction task for central
[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(),
70 fNPhiBins(40)
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(),
93 fNPhiBins(40)
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 //
141 fInspector = o.fInspector;
142 fTrackDensity = o.fTrackDensity;
143 fVtxBins = o.fVtxBins;
144 fFirstEvent = o.fFirstEvent;
145 fHEvents = o.fHEvents;
146 fHEventsTr = o.fHEventsTr;
147 fHEventsTrVtx = o.fHEventsTrVtx;
148 SetVertexAxis(o.fVtxAxis);
149 SetEtaAxis(o.fEtaAxis);
150 fNPhiBins = o.fNPhiBins;
151
152 return *this;
153}
154
155//____________________________________________________________________
156void
157AliCentralMCCorrectionsTask::Init()
158{
159 //
160 // Initialize the task
161 //
162 //
163}
164
165//____________________________________________________________________
166void
167AliCentralMCCorrectionsTask::SetVertexAxis(Int_t nBin, Double_t min,
168 Double_t max)
169{
170 //
171 // Set the vertex axis to use
172 //
173 // Parameters:
174 // nBins Number of bins
175 // vzMin Least @f$z@f$ coordinate of interation point
176 // vzMax Largest @f$z@f$ coordinate of interation point
177 //
178 if (max < min) max = -min;
179 if (min < max) {
180 Double_t tmp = min;
181 min = max;
182 max = tmp;
183 }
184 if (min < -10)
185 AliWarning(Form("Minimum vertex %f < -10, make sure you want this",min));
186 if (max > +10)
187 AliWarning(Form("Minimum vertex %f > +10, make sure you want this",max));
188 fVtxAxis.Set(nBin, min, max);
189}
190//____________________________________________________________________
191void
192AliCentralMCCorrectionsTask::SetVertexAxis(const TAxis& axis)
193{
194 //
195 // Set the vertex axis to use
196 //
197 // Parameters:
198 // axis Axis
199 //
200 SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
201}
202
203//____________________________________________________________________
204void
205AliCentralMCCorrectionsTask::SetEtaAxis(Int_t nBin, Double_t min, Double_t max)
206{
207 //
208 // Set the eta axis to use
209 //
210 // Parameters:
211 // nBins Number of bins
212 // vzMin Least @f$\eta@f$
213 // vzMax Largest @f$\eta@f$
214 //
215 if (max < min) max = -min;
216 if (min < max) {
217 Double_t tmp = min;
218 min = max;
219 max = tmp;
220 }
221 if (min < -4)
222 AliWarning(Form("Minimum eta %f < -4, make sure you want this",min));
223 if (max > +6)
224 AliWarning(Form("Minimum eta %f > +6, make sure you want this",max));
225 fVtxAxis.Set(nBin, min, max);
226}
227//____________________________________________________________________
228void
229AliCentralMCCorrectionsTask::SetEtaAxis(const TAxis& axis)
230{
231 //
232 // Set the eta axis to use
233 //
234 // Parameters:
235 // axis Axis
236 //
237 SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
238}
239
240//____________________________________________________________________
241void
242AliCentralMCCorrectionsTask::DefineBins(TList* l)
243{
244 if (!fVtxBins) fVtxBins = new TObjArray(fVtxAxis.GetNbins(), 1);
245 if (fVtxBins->GetEntries() > 0) return;
246
247 fVtxBins->SetOwner();
248 for (Int_t i = 1; i <= fVtxAxis.GetNbins(); i++) {
249 Double_t low = fVtxAxis.GetBinLowEdge(i);
250 Double_t high = fVtxAxis.GetBinUpEdge(i);
251 VtxBin* bin = new VtxBin(low, high, fEtaAxis, fNPhiBins);
252 fVtxBins->AddAt(bin, i);
253 bin->DefineOutput(l);
254 }
255}
256
257//____________________________________________________________________
258void
259AliCentralMCCorrectionsTask::UserCreateOutputObjects()
260{
261 //
262 // Create output objects
263 //
264 //
265 fList = new TList;
266 fList->SetOwner();
267 fList->SetName(Form("%sSums", GetName()));
268
269 DefineBins(fList);
270
271 fHEvents = new TH1I(GetEventName(false,false),
272 "Number of all events",
273 fVtxAxis.GetNbins(),
274 fVtxAxis.GetXmin(),
275 fVtxAxis.GetXmax());
276 fHEvents->SetXTitle("v_{z} [cm]");
277 fHEvents->SetYTitle("# of events");
278 fHEvents->SetFillColor(kBlue+1);
279 fHEvents->SetFillStyle(3001);
280 fHEvents->SetDirectory(0);
281 fList->Add(fHEvents);
282
283 fHEventsTr = new TH1I(GetEventName(true, false),
284 "Number of triggered events",
285 fVtxAxis.GetNbins(),
286 fVtxAxis.GetXmin(),
287 fVtxAxis.GetXmax());
288 fHEventsTr->SetXTitle("v_{z} [cm]");
289 fHEventsTr->SetYTitle("# of events");
290 fHEventsTr->SetFillColor(kRed+1);
291 fHEventsTr->SetFillStyle(3001);
292 fHEventsTr->SetDirectory(0);
293 fList->Add(fHEventsTr);
294
295 fHEventsTrVtx = new TH1I(GetEventName(true,true),
296 "Number of events w/trigger and vertex",
297 fVtxAxis.GetNbins(),
298 fVtxAxis.GetXmin(),
299 fVtxAxis.GetXmax());
300 fHEventsTrVtx->SetXTitle("v_{z} [cm]");
301 fHEventsTrVtx->SetYTitle("# of events");
302 fHEventsTrVtx->SetFillColor(kBlue+1);
303 fHEventsTrVtx->SetFillStyle(3001);
304 fHEventsTrVtx->SetDirectory(0);
305 fList->Add(fHEventsTrVtx);
306
307 // Copy axis objects to output
308 TH1* vtxAxis = new TH1D("vtxAxis", "Vertex axis",
309 fVtxAxis.GetNbins(),
310 fVtxAxis.GetXmin(),
311 fVtxAxis.GetXmax());
312 TH1* etaAxis = new TH1D("etaAxis", "Eta axis",
313 fEtaAxis.GetNbins(),
314 fEtaAxis.GetXmin(),
315 fEtaAxis.GetXmax());
316 fList->Add(vtxAxis);
317 fList->Add(etaAxis);
318
319 AliInfo(Form("Initialising sub-routines: %p, %p",
320 &fInspector, &fTrackDensity));
321 fInspector.DefineOutput(fList);
322 fInspector.Init(fVtxAxis);
323 fTrackDensity.DefineOutput(fList);
324
325 PostData(1, fList);
326}
327//____________________________________________________________________
328void
329AliCentralMCCorrectionsTask::UserExec(Option_t*)
330{
331 //
332 // Process each event
333 //
334 // Parameters:
335 // option Not used
336 //
337
338 // Get the input data - MC event
339 AliMCEvent* mcEvent = MCEvent();
340 if (!mcEvent) {
341 AliWarning("No MC event found");
342 return;
343 }
344
345 // Get the input data - ESD event
346 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
347 if (!esd) {
348 AliWarning("No ESD event found for input event");
349 return;
350 }
351
352 //--- Read run information -----------------------------------------
353 if (fFirstEvent && esd->GetESDRun()) {
354 fInspector.ReadRunDetails(esd);
355
356 AliInfo(Form("Initializing with parameters from the ESD:\n"
357 " AliESDEvent::GetBeamEnergy() ->%f\n"
358 " AliESDEvent::GetBeamType() ->%s\n"
359 " AliESDEvent::GetCurrentL3() ->%f\n"
360 " AliESDEvent::GetMagneticField()->%f\n"
361 " AliESDEvent::GetRunNumber() ->%d\n",
362 esd->GetBeamEnergy(),
363 esd->GetBeamType(),
364 esd->GetCurrentL3(),
365 esd->GetMagneticField(),
366 esd->GetRunNumber()));
367
368 Print();
369 fFirstEvent = false;
370 }
371
372 // Some variables
373 UInt_t triggers; // Trigger bits
374 Bool_t lowFlux; // Low flux flag
375 UShort_t iVz; // Vertex bin from ESD
376 Double_t vZ; // Z coordinate from ESD
377 Double_t cent; // Centrality
378 UShort_t iVzMc; // Vertex bin from MC
379 Double_t vZMc; // Z coordinate of IP vertex from MC
380 Double_t b; // Impact parameter
381 Int_t nPart; // Number of participants
382 Int_t nBin; // Number of binary collisions
383 Double_t phiR; // Reaction plane from MC
384 UShort_t nClusters;// Number of SPD clusters
385 // Process the data
386 UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ,
387 cent, nClusters);
388 fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc,
389 b, nPart, nBin, phiR);
390
391 Bool_t isInel = triggers & AliAODForwardMult::kInel;
392 Bool_t hasVtx = retESD == AliFMDMCEventInspector::kOk;
393
394 // Fill the event count histograms
395 if (isInel) fHEventsTr->Fill(vZMc);
396 if (isInel && hasVtx) fHEventsTrVtx->Fill(vZMc);
397 fHEvents->Fill(vZMc);
398
399 // Now find our vertex bin object
400 VtxBin* bin = 0;
401 if (iVzMc > 0 && iVzMc <= fVtxAxis.GetNbins())
402 bin = static_cast<VtxBin*>(fVtxBins->At(iVzMc));
403 if (!bin) {
404 // AliError(Form("No vertex bin object @ %d (%f)", iVzMc, vZMc));
405 return;
406 }
407
408 // Now process our input data and store in own ESD object
409 fTrackDensity.Calculate(*mcEvent, vZMc, *bin->fHits, bin->fPrimary);
410 bin->fCounts->Fill(0.5);
411}
412
413//____________________________________________________________________
414void
415AliCentralMCCorrectionsTask::Terminate(Option_t*)
416{
417 //
418 // End of job
419 //
420 // Parameters:
421 // option Not used
422 //
423
424 fList = dynamic_cast<TList*>(GetOutputData(1));
425 if (!fList) {
426 AliError("No output list defined");
427 return;
428 }
429
430 // Output list
431 TList* output = new TList;
432 output->SetOwner();
433 output->SetName(Form("%sResults", GetName()));
434
435 // --- Fill correction object --------------------------------------
436 AliCentralCorrSecondaryMap* corr = new AliCentralCorrSecondaryMap;
437 corr->SetVertexAxis(fVtxAxis);
438 // corr->SetEtaAxis(fEtaAxis);
439
440 TIter next(fVtxBins);
441 VtxBin* bin = 0;
442 UShort_t iVz = 1;
443 while ((bin = static_cast<VtxBin*>(next())))
444 bin->Finish(fList, output, iVz++, corr);
445
446 output->Add(corr);
447
448 PostData(2, output);
449}
450
451//____________________________________________________________________
452void
453AliCentralMCCorrectionsTask::Print(Option_t* option) const
454{
455 std::cout << ClassName() << "\n"
456 << " Vertex bins: " << fVtxAxis.GetNbins() << '\n'
457 << " Vertex range: [" << fVtxAxis.GetXmin()
458 << "," << fVtxAxis.GetXmax() << "]\n"
459 << " Eta bins: " << fEtaAxis.GetNbins() << '\n'
460 << " Eta range: [" << fEtaAxis.GetXmin()
461 << "," << fEtaAxis.GetXmax() << "]\n"
462 << " # of phi bins: " << fNPhiBins
463 << std::endl;
464 gROOT->IncreaseDirLevel();
465 fInspector.Print(option);
466 fTrackDensity.Print(option);
467 gROOT->DecreaseDirLevel();
468}
469
470//====================================================================
471const char*
472AliCentralMCCorrectionsTask::VtxBin::BinName(Double_t low,
473 Double_t high)
474{
475 static TString buf;
476 buf = Form("vtx%+05.1f_%+05.1f", low, high);
477 buf.ReplaceAll("+", "p");
478 buf.ReplaceAll("-", "m");
479 buf.ReplaceAll(".", "d");
480 return buf.Data();
481}
482
483
484//____________________________________________________________________
485AliCentralMCCorrectionsTask::VtxBin::VtxBin()
486 : fHits(0),
487 fPrimary(0),
488 fCounts(0)
489{
490}
491//____________________________________________________________________
492AliCentralMCCorrectionsTask::VtxBin::VtxBin(Double_t low,
493 Double_t high,
494 const TAxis& axis,
495 UShort_t nPhi)
496 : TNamed(BinName(low, high),
497 Form("%+5.1fcm<v_{z}<%+5.1fcm", low, high)),
498 fHits(0),
499 fPrimary(0),
500 fCounts(0)
501{
502 fPrimary = new TH2D("primary", "Primaries",
503 axis.GetNbins(), axis.GetXmin(), axis.GetXmax(),
504 nPhi, 0, 2*TMath::Pi());
505 fPrimary->SetXTitle("#eta");
506 fPrimary->SetYTitle("#varphi [radians]");
507 fPrimary->Sumw2();
508 fPrimary->SetDirectory(0);
509
510 fHits = static_cast<TH2D*>(fPrimary->Clone("hits"));
511 fHits->SetTitle("Hits");
512 fHits->SetDirectory(0);
513
514 fCounts = new TH1D("counts", "Counts", 1, 0, 1);
515 fCounts->SetXTitle("Events");
516 fCounts->SetYTitle("# of Events");
517 fCounts->SetDirectory(0);
518}
519
520//____________________________________________________________________
521AliCentralMCCorrectionsTask::VtxBin::VtxBin(const VtxBin& o)
522 : TNamed(o),
523 fHits(0),
524 fPrimary(0),
525 fCounts(0)
526{
527 if (o.fHits) {
528 fHits = static_cast<TH2D*>(o.fHits->Clone());
529 fHits->SetDirectory(0);
530 }
531 if (o.fPrimary) {
532 fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
533 fPrimary->SetDirectory(0);
534 }
535 if (o.fCounts) {
536 fCounts = static_cast<TH1D*>(o.fCounts->Clone());
537 fCounts->SetDirectory(0);
538 }
539}
540
541//____________________________________________________________________
542AliCentralMCCorrectionsTask::VtxBin&
543AliCentralMCCorrectionsTask::VtxBin::operator=(const VtxBin& o)
544{
545 TNamed::operator=(o);
546 fHits = 0;
547 fPrimary = 0;
548 fCounts = 0;
549 if (o.fHits) {
550 fHits = static_cast<TH2D*>(o.fHits->Clone());
551 fHits->SetDirectory(0);
552 }
553 if (o.fPrimary) {
554 fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
555 fPrimary->SetDirectory(0);
556 }
557 if (o.fCounts) {
558 fCounts = static_cast<TH1D*>(o.fCounts->Clone());
559 fCounts->SetDirectory(0);
560 }
561 return *this;
562}
563
564//____________________________________________________________________
565void
566AliCentralMCCorrectionsTask::VtxBin::DefineOutput(TList* l)
567{
568 TList* d = new TList;
569 d->SetName(GetName());
570 d->SetOwner();
571 l->Add(d);
572
573 d->Add(fHits);
574 d->Add(fPrimary);
575 d->Add(fCounts);
576}
577//____________________________________________________________________
578void
579AliCentralMCCorrectionsTask::VtxBin::Finish(const TList* input,
580 TList* output,
581 UShort_t iVz,
582 AliCentralCorrSecondaryMap* map)
583{
584 TList* out = new TList;
585 out->SetName(GetName());
586 out->SetOwner();
587 output->Add(out);
588
589 TList* l = static_cast<TList*>(input->FindObject(GetName()));
590 if (!l) {
591 AliError(Form("List %s not found in %s", GetName(), input->GetName()));
592 return;
593 }
594
595
596 TH2D* hits = static_cast<TH2D*>(l->FindObject("hits"));
597 TH2D* prim = static_cast<TH2D*>(l->FindObject("primary"));
598 if (!hits || !prim) {
599 AliError(Form("Missing histograms: %p, %p", hits, prim));
600 return;
601 }
602
603 TH2D* h = static_cast<TH2D*>(hits->Clone("bgCorr"));
604 h->SetDirectory(0);
605 h->Divide(fPrimary);
606
607 map->SetCorrection(iVz, h);
608
609 out->Add(h);
610}
611
612//
613// EOF
614//