2 // Calculate the corrections in the forward regions
14 #include "AliForwardMCCorrectionsTask.h"
15 #include "AliTriggerAnalysis.h"
16 #include "AliPhysicsSelection.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"
25 #include "AliMCEvent.h"
26 #include "AliAODForwardMult.h"
27 #include "AliFMDStripIndex.h"
28 #include "AliFMDCorrSecondaryMap.h"
31 #include <TDirectory.h>
37 //====================================================================
39 const char* GetEventName(Bool_t tr, Bool_t vtx)
41 return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : ""));
44 const char* GetHitsName(UShort_t d, Char_t r)
46 return Form("hitsFMD%d%c", d, r);
48 const char* GetStripsName(UShort_t d, Char_t r)
50 return Form("stripsFMD%d%c", d, r);
52 const char* GetPrimaryName(Char_t r, Bool_t trVtx)
54 return Form("primaries%s%s",
55 ((r == 'I' || r == 'i') ? "Inner" : "Outer"),
56 (trVtx ? "TrVtx" : "All"));
61 //====================================================================
62 AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask()
63 : AliAnalysisTaskSE(),
84 //____________________________________________________________________
85 AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const char* name)
86 : AliAnalysisTaskSE(name),
87 fInspector("eventInspector"),
88 fTrackDensity("trackDensity"),
105 DefineOutput(1, TList::Class());
106 DefineOutput(2, TList::Class());
108 "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
109 "AliESDFMD.,SPDVertex.,PrimaryVertex.";
112 //____________________________________________________________________
113 AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const AliForwardMCCorrectionsTask& o)
114 : AliAnalysisTaskSE(o),
115 fInspector(o.fInspector),
119 fFirstEvent(o.fFirstEvent),
120 fHEvents(o.fHEvents),
121 fHEventsTr(o.fHEventsTr),
122 fHEventsTrVtx(o.fHEventsTrVtx),
131 // o Object to copy from
135 //____________________________________________________________________
136 AliForwardMCCorrectionsTask&
137 AliForwardMCCorrectionsTask::operator=(const AliForwardMCCorrectionsTask& o)
140 // Assignment operator
143 // o Object to assign from
146 // Reference to this object
148 if (&o == this) return *this;
149 fInspector = o.fInspector;
150 fTrackDensity = o.fTrackDensity;
152 fVtxBins = o.fVtxBins;
153 fFirstEvent = o.fFirstEvent;
154 fHEvents = o.fHEvents;
155 fHEventsTr = o.fHEventsTr;
156 fHEventsTrVtx = o.fHEventsTrVtx;
157 SetVertexAxis(o.fVtxAxis);
158 SetEtaAxis(o.fEtaAxis);
163 //____________________________________________________________________
165 AliForwardMCCorrectionsTask::Init()
168 // Initialize the task
173 //____________________________________________________________________
175 AliForwardMCCorrectionsTask::SetVertexAxis(Int_t nBin, Double_t min,
179 // Set the vertex axis to use
182 // nBins Number of bins
183 // vzMin Least @f$z@f$ coordinate of interation point
184 // vzMax Largest @f$z@f$ coordinate of interation point
193 AliWarning(Form("Minimum vertex %f < -10, make sure you want this",min));
195 AliWarning(Form("Minimum vertex %f > +10, make sure you want this",max));
197 fVtxAxis.Set(nBin, min, max);
199 //____________________________________________________________________
201 AliForwardMCCorrectionsTask::SetVertexAxis(const TAxis& axis)
204 // Set the vertex axis to use
209 SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
212 //____________________________________________________________________
214 AliForwardMCCorrectionsTask::SetEtaAxis(Int_t nBin, Double_t min, Double_t max)
217 // Set the eta axis to use
220 // nBins Number of bins
221 // vzMin Least @f$\eta@f$
222 // vzMax Largest @f$\eta@f$
230 AliWarning(Form("Minimum eta %f < -4, make sure you want this",min));
232 AliWarning(Form("Minimum eta %f > +6, make sure you want this",max));
233 fEtaAxis.Set(nBin, min, max);
235 //____________________________________________________________________
237 AliForwardMCCorrectionsTask::SetEtaAxis(const TAxis& axis)
240 // Set the eta axis to use
245 SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
248 //____________________________________________________________________
250 AliForwardMCCorrectionsTask::DefineBins(TList* l)
252 if (!fVtxBins) fVtxBins = new TObjArray(fVtxAxis.GetNbins(), 1);
253 if (fVtxBins->GetEntries() > 0) return;
255 fVtxBins->SetOwner();
256 for (Int_t i = 1; i <= fVtxAxis.GetNbins(); i++) {
257 Double_t low = fVtxAxis.GetBinLowEdge(i);
258 Double_t high = fVtxAxis.GetBinUpEdge(i);
259 VtxBin* bin = new VtxBin(low, high, fEtaAxis);
260 fVtxBins->AddAt(bin, i);
261 bin->DefineOutput(l);
265 //____________________________________________________________________
267 AliForwardMCCorrectionsTask::UserCreateOutputObjects()
270 // Create output objects
275 fList->SetName(Form("%sSums", GetName()));
279 fHEvents = new TH1I(GetEventName(false,false),
280 "Number of all events",
284 fHEvents->SetXTitle("v_{z} [cm]");
285 fHEvents->SetYTitle("# of events");
286 fHEvents->SetFillColor(kBlue+1);
287 fHEvents->SetFillStyle(3001);
288 fHEvents->SetDirectory(0);
289 fList->Add(fHEvents);
291 fHEventsTr = new TH1I(GetEventName(true, false),
292 "Number of triggered events",
296 fHEventsTr->SetXTitle("v_{z} [cm]");
297 fHEventsTr->SetYTitle("# of events");
298 fHEventsTr->SetFillColor(kRed+1);
299 fHEventsTr->SetFillStyle(3001);
300 fHEventsTr->SetDirectory(0);
301 fList->Add(fHEventsTr);
303 fHEventsTrVtx = new TH1I(GetEventName(true,true),
304 "Number of events w/trigger and vertex",
308 fHEventsTrVtx->SetXTitle("v_{z} [cm]");
309 fHEventsTrVtx->SetYTitle("# of events");
310 fHEventsTrVtx->SetFillColor(kBlue+1);
311 fHEventsTrVtx->SetFillStyle(3001);
312 fHEventsTrVtx->SetDirectory(0);
313 fList->Add(fHEventsTrVtx);
315 // Copy axis objects to output
316 TH1* vtxAxis = new TH1D("vtxAxis", "Vertex axis",
320 TH1* etaAxis = new TH1D("etaAxis", "Eta axis",
327 AliInfo(Form("Initialising sub-routines: %p, %p",
328 &fInspector, &fTrackDensity));
329 fInspector.DefineOutput(fList);
330 fTrackDensity.DefineOutput(fList);
334 //____________________________________________________________________
336 AliForwardMCCorrectionsTask::UserExec(Option_t*)
339 // Process each event
345 // Get the input data - MC event
346 AliMCEvent* mcEvent = MCEvent();
348 AliWarning("No MC event found");
352 // Get the input data - ESD event
353 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
355 AliWarning("No ESD event found for input event");
359 //--- Read run information -----------------------------------------
360 if (fFirstEvent && esd->GetESDRun()) {
361 fInspector.ReadRunDetails(esd);
363 AliInfo(Form("Initializing with parameters from the ESD:\n"
364 " AliESDEvent::GetBeamEnergy() ->%f\n"
365 " AliESDEvent::GetBeamType() ->%s\n"
366 " AliESDEvent::GetCurrentL3() ->%f\n"
367 " AliESDEvent::GetMagneticField()->%f\n"
368 " AliESDEvent::GetRunNumber() ->%d\n",
369 esd->GetBeamEnergy(),
372 esd->GetMagneticField(),
373 esd->GetRunNumber()));
375 fInspector.Init(fVtxAxis);
382 UInt_t triggers; // Trigger bits
383 Bool_t lowFlux; // Low flux flag
384 UShort_t iVz; // Vertex bin from ESD
385 Double_t vZ; // Z coordinate from ESD
386 Double_t cent; // Centrality
387 UShort_t iVzMc; // Vertex bin from MC
388 Double_t vZMc; // Z coordinate of IP vertex from MC
389 Double_t b; // Impact parameter
390 Double_t cMC; // Centrality estimate from b
391 Int_t nPart; // Number of participants
392 Int_t nBin; // Number of binary collisions
393 Double_t phiR; // Reaction plane from MC
394 UShort_t nClusters;// Number of SPD clusters
396 UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ,
398 fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc,
399 b, cMC, nPart, nBin, phiR);
401 Bool_t isInel = triggers & AliAODForwardMult::kInel;
402 Bool_t hasVtx = retESD == AliFMDMCEventInspector::kOk;
404 // Fill the event count histograms
405 if (isInel) fHEventsTr->Fill(vZMc);
406 if (isInel && hasVtx) fHEventsTrVtx->Fill(vZMc);
407 fHEvents->Fill(vZMc);
409 // Now find our vertex bin object
411 if (iVzMc > 0 && iVzMc <= fVtxAxis.GetNbins())
412 bin = static_cast<VtxBin*>(fVtxBins->At(iVzMc));
414 // AliError(Form("No vertex bin object @ %d (%f)", iVzMc, vZMc));
418 // Clear our ESD object
422 AliESDFMD* esdFMD = esd->GetFMDData();
424 // Now process our input data and store in own ESD object
425 fTrackDensity.Calculate(*esdFMD, *mcEvent, vZMc, fESDFMD, bin->fPrimary);
426 bin->fCounts->Fill(0.5);
428 // And then bin the data in our vtxbin
429 for (UShort_t d=1; d<=3; d++) {
430 UShort_t nr = (d == 1 ? 1 : 2);
431 for (UShort_t q=0; q<nr; q++) {
432 Char_t r = (q == 0 ? 'I' : 'O');
433 UShort_t ns= (q == 0 ? 20 : 40);
434 UShort_t nt= (q == 0 ? 512 : 256);
435 TH2D* h = bin->fHists.Get(d,r);
437 for (UShort_t s=0; s<ns; s++) {
438 for (UShort_t t=0; t<nt; t++) {
439 Float_t mult = fESDFMD.Multiplicity(d,r,s,t);
441 if (mult == 0 || mult > 20) continue;
443 Float_t phi = fESDFMD.Phi(d,r,s,t) / 180 * TMath::Pi();
444 Float_t eta = fESDFMD.Eta(d,r,s,t);
445 h->Fill(eta,phi,mult);
453 //____________________________________________________________________
455 AliForwardMCCorrectionsTask::Terminate(Option_t*)
464 fList = dynamic_cast<TList*>(GetOutputData(1));
466 AliError("No output list defined");
473 TList* output = new TList;
475 output->SetName(Form("%sResults", GetName()));
477 // --- Fill correction object --------------------------------------
478 AliFMDCorrSecondaryMap* corr = new AliFMDCorrSecondaryMap;
479 corr->SetVertexAxis(fVtxAxis);
480 corr->SetEtaAxis(fEtaAxis);
482 TIter next(fVtxBins);
485 while ((bin = static_cast<VtxBin*>(next())))
486 bin->Finish(fList, output, iVz++, corr);
493 //____________________________________________________________________
495 AliForwardMCCorrectionsTask::Print(Option_t* option) const
497 std::cout << ClassName() << "\n"
498 << " Vertex bins: " << fVtxAxis.GetNbins() << '\n'
499 << " Vertex range: [" << fVtxAxis.GetXmin()
500 << "," << fVtxAxis.GetXmax() << "]\n"
501 << " Eta bins: " << fEtaAxis.GetNbins() << '\n'
502 << " Eta range: [" << fEtaAxis.GetXmin()
503 << "," << fEtaAxis.GetXmax() << "]"
505 gROOT->IncreaseDirLevel();
506 fInspector.Print(option);
507 fTrackDensity.Print(option);
508 gROOT->DecreaseDirLevel();
511 //====================================================================
513 AliForwardMCCorrectionsTask::VtxBin::BinName(Double_t low,
517 buf = Form("vtx%+05.1f_%+05.1f", low, high);
518 buf.ReplaceAll("+", "p");
519 buf.ReplaceAll("-", "m");
520 buf.ReplaceAll(".", "d");
525 //____________________________________________________________________
526 AliForwardMCCorrectionsTask::VtxBin::VtxBin()
532 //____________________________________________________________________
533 AliForwardMCCorrectionsTask::VtxBin::VtxBin(Double_t low,
536 : TNamed(BinName(low, high),
537 Form("%+5.1fcm<v_{z}<%+5.1fcm", low, high)),
544 fPrimary = new TH2D("primary", "Primaries",
545 axis.GetNbins(), axis.GetXmin(), axis.GetXmax(),
546 40, 0, 2*TMath::Pi());
547 fPrimary->SetXTitle("#eta");
548 fPrimary->SetYTitle("#varphi [radians]");
550 fPrimary->SetDirectory(0);
552 fCounts = new TH1D("counts", "Counts", 1, 0, 1);
553 fCounts->SetXTitle("Events");
554 fCounts->SetYTitle("# of Events");
555 fCounts->SetDirectory(0);
558 //____________________________________________________________________
559 AliForwardMCCorrectionsTask::VtxBin::VtxBin(const VtxBin& o)
566 fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
567 fPrimary->SetDirectory(0);
570 fCounts = static_cast<TH1D*>(o.fCounts->Clone());
571 fCounts->SetDirectory(0);
575 //____________________________________________________________________
576 AliForwardMCCorrectionsTask::VtxBin&
577 AliForwardMCCorrectionsTask::VtxBin::operator=(const VtxBin& o)
579 if (&o == this) return *this;
580 TNamed::operator=(o);
585 fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
586 fPrimary->SetDirectory(0);
589 fCounts = static_cast<TH1D*>(o.fCounts->Clone());
590 fCounts->SetDirectory(0);
595 //____________________________________________________________________
597 AliForwardMCCorrectionsTask::VtxBin::DefineOutput(TList* l)
599 TList* d = new TList;
600 d->SetName(GetName());
604 d->Add(fHists.fFMD1i);
605 d->Add(fHists.fFMD2i);
606 d->Add(fHists.fFMD2o);
607 d->Add(fHists.fFMD3i);
608 d->Add(fHists.fFMD3o);
614 //____________________________________________________________________
616 AliForwardMCCorrectionsTask::VtxBin::MakeBg(const TH2D* hits,
617 const TH2D* primary) const
619 TH2D* h = static_cast<TH2D*>(hits->Clone());
621 TString n(h->GetName());
622 n.ReplaceAll("_cache", "");
629 //____________________________________________________________________
631 AliForwardMCCorrectionsTask::VtxBin::Finish(const TList* input,
634 AliFMDCorrSecondaryMap* map)
636 TList* out = new TList;
637 out->SetName(GetName());
641 TList* l = static_cast<TList*>(input->FindObject(GetName()));
643 AliError(Form("List %s not found in %s", GetName(), input->GetName()));
647 TH2D* fmd1i = static_cast<TH2D*>(l->FindObject("FMD1I_cache"));
648 TH2D* fmd2i = static_cast<TH2D*>(l->FindObject("FMD2I_cache"));
649 TH2D* fmd2o = static_cast<TH2D*>(l->FindObject("FMD2O_cache"));
650 TH2D* fmd3i = static_cast<TH2D*>(l->FindObject("FMD3I_cache"));
651 TH2D* fmd3o = static_cast<TH2D*>(l->FindObject("FMD3O_cache"));
652 TH2D* primO = static_cast<TH2D*>(l->FindObject("primary"));
653 if (!fmd1i || !fmd2i || !fmd2o || !fmd3i || !fmd3o || !primO) {
654 AliError(Form("Missing histogram(s): %p,%p,%p,%p,%p,%p",
655 fmd1i, fmd2i, fmd2o, fmd3i, fmd3o, primO));
659 // Half coverage in phi for inners
660 TH2D* primI = static_cast<TH2D*>(primO->Clone());
661 primI->SetDirectory(0);
664 TH2D* bg1i = MakeBg(fmd1i, primI);
665 TH2D* bg2i = MakeBg(fmd2i, primI);
666 TH2D* bg2o = MakeBg(fmd2o, primO);
667 TH2D* bg3i = MakeBg(fmd3i, primI);
668 TH2D* bg3o = MakeBg(fmd3o, primO);
669 map->SetCorrection(1, 'I', iVz, bg1i);
670 map->SetCorrection(2, 'I', iVz, bg2i);
671 map->SetCorrection(2, 'O', iVz, bg2o);
672 map->SetCorrection(3, 'I', iVz, bg3i);
673 map->SetCorrection(3, 'O', iVz, bg3o);