1 //====================================================================
2 #include "AliBasedNdetaTask.h"
8 #include <TParameter.h>
9 #include <AliAnalysisManager.h>
10 #include <AliAODEvent.h>
11 #include <AliAODHandler.h>
12 #include <AliAODInputHandler.h>
13 #include "AliForwardUtil.h"
14 #include "AliAODForwardMult.h"
18 //____________________________________________________________________
19 AliBasedNdetaTask::AliBasedNdetaTask()
20 : AliAnalysisTaskSE(),
21 fSums(0), // Container of sums
22 fOutput(0), // Container of output
23 fVtxMin(0), // Minimum v_z
24 fVtxMax(0), // Maximum v_z
25 fTriggerMask(0), // Trigger mask
26 fRebin(0), // Rebinning factor
33 fListOfCentralities(0),
38 fNormalizationScheme(kFull),
46 DGUARD(0,0,"Default construction of AliBasedNdetaTask");
49 //____________________________________________________________________
50 AliBasedNdetaTask::AliBasedNdetaTask(const char* name)
51 : AliAnalysisTaskSE(name),
52 fSums(0), // Container of sums
53 fOutput(0), // Container of output
54 fVtxMin(-10), // Minimum v_z
55 fVtxMax(10), // Maximum v_z
56 fTriggerMask(AliAODForwardMult::kInel),
57 fRebin(5), // Rebinning factor
64 fListOfCentralities(0),
69 fNormalizationScheme(kFull),
77 DGUARD(0,0,"Named construction of AliBasedNdetaTask: %s", name);
78 fListOfCentralities = new TObjArray(1);
80 // Set the normalisation scheme
81 SetNormalizationScheme(kFull);
83 // Set the trigger mask
84 SetTriggerMask(AliAODForwardMult::kInel);
86 // Output slot #1 writes into a TH1 container
87 DefineOutput(1, TList::Class());
88 DefineOutput(2, TList::Class());
91 //____________________________________________________________________
92 AliBasedNdetaTask::AliBasedNdetaTask(const AliBasedNdetaTask& o)
93 : AliAnalysisTaskSE(o),
94 fSums(o.fSums), // TList* - Container of sums
95 fOutput(o.fOutput), // Container of output
96 fVtxMin(o.fVtxMin), // Double_t - Minimum v_z
97 fVtxMax(o.fVtxMax), // Double_t - Maximum v_z
98 fTriggerMask(o.fTriggerMask),// Int_t - Trigger mask
99 fRebin(o.fRebin), // Int_t - Rebinning factor
100 fCutEdges(o.fCutEdges), // Bool_t - Whether to cut edges when rebinning
101 fSymmetrice(o.fSymmetrice),
102 fCorrEmpty(o.fCorrEmpty),
103 fUseROOTProj(o.fUseROOTProj),
104 fTriggerEff(o.fTriggerEff),
105 fShapeCorr(o.fShapeCorr),
106 fListOfCentralities(o.fListOfCentralities),
107 fSNNString(o.fSNNString),
108 fSysString(o.fSysString),
110 fCentAxis(o.fCentAxis),
111 fNormalizationScheme(o.fNormalizationScheme),
112 fSchemeString(o.fSchemeString),
113 fTriggerString(o.fTriggerString),
114 fFinalMCCorrFile(o.fFinalMCCorrFile)
116 DGUARD(0,0,"Copy construction of AliBasedNdetaTask");
119 //____________________________________________________________________
120 AliBasedNdetaTask::~AliBasedNdetaTask()
125 DGUARD(fDebug,3,"Destruction of AliBasedNdetaTask");
138 //________________________________________________________________________
140 AliBasedNdetaTask::SetDebugLevel(Int_t lvl)
142 AliAnalysisTaskSE::SetDebugLevel(lvl);
143 for (Int_t i = 0; i < fListOfCentralities->GetEntries(); i++) {
145 static_cast<CentralityBin*>(fListOfCentralities->At(i));
146 bin->SetDebugLevel(lvl);
150 //________________________________________________________________________
152 AliBasedNdetaTask::SetCentralityAxis(UShort_t n, Short_t* bins)
154 DGUARD(fDebug,3,"Set centrality axis, %d bins", n);
156 fCentAxis = new TAxis();
157 fCentAxis->SetName("centAxis");
158 fCentAxis->SetTitle("Centrality [%]");
161 for (UShort_t i = 0; i <= n; i++)
162 dbins[i] = (bins[i] == 100 ? 100.1 : bins[i]);
163 fCentAxis->Set(n, dbins.GetArray());
166 //________________________________________________________________________
168 AliBasedNdetaTask::AddCentralityBin(UShort_t at, Short_t low, Short_t high)
171 // Add a centrality bin
177 DGUARD(fDebug,3,"Add a centrality bin [%f,%f] @ %d", low, high, at);
178 CentralityBin* bin = MakeCentralityBin(GetName(), low, high);
179 bin->SetDebugLevel(fDebug);
180 fListOfCentralities->AddAtAndExpand(bin, at);
183 //________________________________________________________________________
184 AliBasedNdetaTask::CentralityBin*
185 AliBasedNdetaTask::MakeCentralityBin(const char* name,
186 Short_t low, Short_t high) const
189 // Make a centrality bin
192 // name Name used for histograms
193 // low Low cut in percent
194 // high High cut in percent
197 // A newly created centrality bin
199 DGUARD(fDebug,3,"Make a centrality bin [%f,%f]: %s", low, high, name);
200 return new CentralityBin(name, low, high);
203 #define TESTAPPEND(SCHEME,BIT,STRING) \
204 do { if (!(SCHEME & BIT)) break; \
205 if (!s.IsNull()) s.Append(","); s.Append(STRING); } while(false)
207 //________________________________________________________________________
209 AliBasedNdetaTask::NormalizationSchemeString(UShort_t scheme)
211 // Create a string from normalization scheme bits
217 if (scheme == kFull) {
221 TESTAPPEND(scheme, kEventLevel, "EVENT");
222 TESTAPPEND(scheme, kShape, "SHAPE");
223 TESTAPPEND(scheme, kBackground, "BACKGROUND");
224 TESTAPPEND(scheme, kTriggerEfficiency, "TRIGGER");
225 TESTAPPEND(scheme, kZeroBin, "ZEROBIN");
229 //________________________________________________________________________
231 AliBasedNdetaTask::ParseNormalizationScheme(const char* what)
237 TIter next(twhat.Tokenize(" ,|"));
238 while ((opt = static_cast<TObjString*>(next()))) {
239 TString s(opt->GetString());
240 if (s.IsNull()) continue;
243 case '-': add = false; // Fall through
244 case '+': s.Remove(0,1); // Remove character
247 if (s.CompareTo("EVENT") == 0) bit = kEventLevel;
248 else if (s.CompareTo("SHAPE") == 0) bit = kShape;
249 else if (s.CompareTo("BACKGROUND")== 0) bit = kBackground;
250 else if (s.CompareTo("TRIGGER") == 0) bit = kTriggerEfficiency;
251 else if (s.CompareTo("FULL") == 0) bit = kFull;
252 else if (s.CompareTo("NONE") == 0) bit = kNone;
253 else if (s.CompareTo("ZEROBIN") == 0) bit = kZeroBin;
255 ::Warning("SetNormalizationScheme", "Unknown option %s", s.Data());
256 if (add) scheme |= bit;
261 //________________________________________________________________________
263 AliBasedNdetaTask::SetNormalizationScheme(const char* what)
266 // Set normalisation scheme
268 DGUARD(fDebug,3,"Set the normalization scheme: %s", what);
269 SetNormalizationScheme(ParseNormalizationScheme(what));
271 //________________________________________________________________________
273 AliBasedNdetaTask::SetNormalizationScheme(UShort_t scheme)
275 DGUARD(fDebug,3,"Set the normalization scheme: 0x%x", scheme);
276 fNormalizationScheme = scheme;
278 const Char_t* tit = NormalizationSchemeString(scheme);
279 if (!fSchemeString) fSchemeString = new TNamed("scheme", "");
280 fSchemeString->SetTitle(tit);
281 fSchemeString->SetUniqueID(fNormalizationScheme);
284 fSchemeString = new TParameter<int>("scheme", scheme);
285 fSchemeString->SetUniqueID(fNormalizationScheme);
288 //________________________________________________________________________
290 AliBasedNdetaTask::SetTriggerMask(const char* mask)
293 // Set the trigger maskl
298 DGUARD(fDebug,3,"Set the trigger mask: %s", mask);
299 SetTriggerMask(AliAODForwardMult::MakeTriggerMask(mask));
301 //________________________________________________________________________
303 AliBasedNdetaTask::SetTriggerMask(UShort_t mask)
305 DGUARD(fDebug,3,"Set the trigger mask: 0x%0x", mask);
308 TString tit(AliAODForwardMult::GetTriggerString(mask));
309 tit = tit.Strip(TString::kBoth);
310 if (!fTriggerString) fTriggerString = new TNamed("trigger", "");
311 fTriggerString->SetTitle(tit);
312 fTriggerString->SetUniqueID(fTriggerMask);
315 fTriggerString = new TParameter<int>("trigger", fTriggerMask);
316 fTriggerString->SetUniqueID(fTriggerMask);
320 //________________________________________________________________________
322 AliBasedNdetaTask::SetShapeCorrection(const TH2F* c)
325 // Set the shape correction (a.k.a., track correction) for selected
331 DGUARD(fDebug,3,"Set the shape correction: %p", c);
333 fShapeCorr = static_cast<TH2F*>(c->Clone());
334 fShapeCorr->SetDirectory(0);
337 //________________________________________________________________________
339 AliBasedNdetaTask::InitializeCentBins()
341 if (fListOfCentralities->GetEntries() > 0) return;
343 // Automatically add 'all' centrality bin if nothing has been defined.
344 AddCentralityBin(0, 0, 0);
345 if (fCentAxis && fCentAxis->GetNbins() > 0 && fCentAxis->GetXbins()) {
346 const TArrayD* bins = fCentAxis->GetXbins();
347 Int_t nbin = fCentAxis->GetNbins();
348 for (Int_t i = 0; i < nbin; i++)
349 AddCentralityBin(i+1, Short_t((*bins)[i]), Short_t((*bins)[i+1]));
353 //________________________________________________________________________
355 AliBasedNdetaTask::UserCreateOutputObjects()
358 // Create output objects.
360 // This is called once per slave process
362 DGUARD(fDebug,1,"Create user ouput object");
364 fSums->SetName(Form("%s_sums", GetName()));
367 InitializeCentBins();
368 if (fCentAxis) fSums->Add(fCentAxis);
371 // Centrality histogram
372 fCent = new TH1D("cent", "Centrality", 100, 0, 100);
373 fCent->SetDirectory(0);
377 // Loop over centrality bins
378 TIter next(fListOfCentralities);
379 CentralityBin* bin = 0;
380 while ((bin = static_cast<CentralityBin*>(next())))
381 bin->CreateOutputObjects(fSums);
384 // Post data for ALL output slots >0 here, to get at least an empty
389 //____________________________________________________________________
391 AliBasedNdetaTask::UserExec(Option_t *)
394 // Process a single event
400 DGUARD(fDebug,1,"Analyse the AOD event");
401 AliAODEvent* aod = AliForwardUtil::GetAODEvent(this);
405 TObject* obj = aod->FindListObject("Forward");
407 AliWarning("No forward object found");
410 AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
412 // Fill centrality histogram
413 Float_t cent = forward->GetCentrality();
416 // Get the histogram(s)
417 TH2D* data = GetHistogram(aod, false);
418 TH2D* dataMC = GetHistogram(aod, true);
420 Bool_t isZero = ((fNormalizationScheme & kZeroBin) &&
421 !forward->IsTriggerBits(AliAODForwardMult::kNClusterGt0));
424 // Loop over centrality bins
425 CentralityBin* allBin =
426 static_cast<CentralityBin*>(fListOfCentralities->At(0));
427 allBin->ProcessEvent(forward, fTriggerMask, isZero,
428 fVtxMin, fVtxMax, data, dataMC);
430 // Find this centrality bin
431 if (fCentAxis && fCentAxis->GetNbins() > 0) {
432 Int_t icent = fCentAxis->FindBin(cent);
433 CentralityBin* thisBin = 0;
434 if (icent >= 1 && icent <= fCentAxis->GetNbins())
435 thisBin = static_cast<CentralityBin*>(fListOfCentralities->At(icent));
437 thisBin->ProcessEvent(forward, fTriggerMask, isZero, fVtxMin, fVtxMax,
441 // Here, we get the update
444 UShort_t sNN = forward->GetSNN();
445 fSNNString = new TNamed("sNN", "");
446 fSNNString->SetTitle(AliForwardUtil::CenterOfMassEnergyString(sNN));
447 fSNNString->SetUniqueID(sNN);
448 fSums->Add(fSNNString);
450 UShort_t sys = forward->GetSystem();
451 fSysString = new TNamed("sys", "");
452 fSysString->SetTitle(AliForwardUtil::CollisionSystemString(sys));
453 fSysString->SetUniqueID(sys);
454 fSums->Add(fSysString);
456 fSums->Add(fSchemeString);
457 fSums->Add(fTriggerString);
463 UShort_t sNN = forward->GetSNN();
464 fSNNString = new TParameter<int>("sNN", sNN);
465 fSNNString->SetUniqueID(sNN);
466 fSums->Add(fSNNString);
468 UShort_t sys = forward->GetSystem();
469 fSysString = new TParameter<int>("sys", sys);
470 fSysString->SetUniqueID(sys);
471 fSums->Add(fSysString);
473 fSums->Add(fSchemeString);
474 fSums->Add(fTriggerString);
483 //________________________________________________________________________
485 AliBasedNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker,
486 const char* title, const char* ytitle)
489 // Set histogram graphical options, etc.
492 // h Histogram to modify
493 // colour Marker color
494 // marker Marker style
495 // title Title of histogram
496 // ytitle Title on y-axis.
499 h->SetMarkerColor(colour);
500 h->SetMarkerStyle(marker);
501 h->SetMarkerSize(marker == 29 || marker == 30 ? 1.2 : 1);
503 h->SetYTitle(ytitle);
504 h->GetXaxis()->SetTitleFont(132);
505 h->GetXaxis()->SetLabelFont(132);
506 h->GetXaxis()->SetNdivisions(10);
507 h->GetYaxis()->SetTitleFont(132);
508 h->GetYaxis()->SetLabelFont(132);
509 h->GetYaxis()->SetNdivisions(10);
510 h->GetYaxis()->SetDecimals();
514 //________________________________________________________________________
516 AliBasedNdetaTask::ScaleToCoverage(TH2D* copy, const TH1D* norm)
518 // Normalize to the acceptance -
519 // dndeta->Divide(accNorm);
520 for (Int_t i = 1; i <= copy->GetNbinsX(); i++) {
521 Double_t a = norm->GetBinContent(i);
522 for (Int_t j = 1; j <= copy->GetNbinsY(); j++) {
524 copy->SetBinContent(i,j,0);
525 copy->SetBinError(i,j,0);
528 Double_t c = copy->GetBinContent(i, j);
529 Double_t e = copy->GetBinError(i, j);
530 copy->SetBinContent(i, j, c / a);
531 copy->SetBinError(i, j, e / a);
536 //________________________________________________________________________
538 AliBasedNdetaTask::ProjectX(const TH2D* h,
547 // Project onto the X axis
552 // firstbin First bin to use
553 // lastbin Last bin to use
554 // error Whether to calculate errors
557 // Newly created histogram or null
561 return h->ProjectionX(name, firstbin, lastbin, (error ? "e" : ""));
563 TAxis* xaxis = h->GetXaxis();
564 TAxis* yaxis = h->GetYaxis();
565 TH1D* ret = new TH1D(name, h->GetTitle(), xaxis->GetNbins(),
566 xaxis->GetXmin(), xaxis->GetXmax());
567 static_cast<const TAttLine*>(h)->Copy(*ret);
568 static_cast<const TAttFill*>(h)->Copy(*ret);
569 static_cast<const TAttMarker*>(h)->Copy(*ret);
570 ret->GetXaxis()->ImportAttributes(xaxis);
572 Int_t first = firstbin;
573 Int_t last = lastbin;
574 if (first < 0) first = 0;
575 else if (first >= yaxis->GetNbins()+1) first = yaxis->GetNbins();
576 if (last < 0) last = yaxis->GetNbins();
577 else if (last > yaxis->GetNbins()+1) last = yaxis->GetNbins();
578 if (last-first < 0) {
579 AliWarningGeneral("AliBasedNdetaTask",
580 Form("Nothing to project [%d,%d]", first, last));
586 // AliInfo(Form("Projecting bins [%d,%d] of %s", first, last, h->GetName()));
587 Int_t ybins = (last-first+1);
588 for (Int_t xbin = 0; xbin <= xaxis->GetNbins()+1; xbin++) {
589 Double_t content = 0;
594 for (Int_t ybin = first; ybin <= last; ybin++) {
595 Double_t c1 = h->GetCellContent(xbin, ybin);
596 Double_t e1 = h->GetCellError(xbin, ybin);
599 if (c1 < 1e-12) continue;
609 if(content > 0 && nbins > 0) {
610 Double_t factor = (corr ? Double_t(ybins) / nbins : 1);
612 AliWarningGeneral(ret->GetName(),
613 Form("factor @ %d is %d/%d -> %f",
614 xbin, ybins, nbins, factor));
617 // Calculate weighted average
618 ret->SetBinContent(xbin, content * factor);
619 ret->SetBinError(xbin, factor * TMath::Sqrt(error2));
622 ret->SetBinContent(xbin, factor * content);
629 //________________________________________________________________________
631 AliBasedNdetaTask::Terminate(Option_t *)
634 // Called at end of event processing..
636 // This is called once in the master
641 // Draw result to screen, or perform fitting, normalizations Called
642 // once at the end of the query
643 DGUARD(fDebug,1,"Process final merged results");
645 fSums = dynamic_cast<TList*> (GetOutputData(1));
647 AliError("Could not retrieve TList fSums");
652 fOutput->SetName(Form("%s_result", GetName()));
655 fSNNString = static_cast<TParameter<int>*>(fSums->FindObject("sNN"));
656 fSysString = static_cast<TParameter<int>*>(fSums->FindObject("sys"));
657 fCentAxis = static_cast<TAxis*>(fSums->FindObject("centAxis"));
658 fSchemeString = static_cast<TParameter<int>*>(fSums->FindObject("scheme"));
659 fTriggerString = static_cast<TParameter<int>*>(fSums->FindObject("trigger"));
661 if(fSysString && fSNNString &&
662 fSysString->GetUniqueID() == AliForwardUtil::kPP)
663 LoadNormalizationData(fSysString->GetUniqueID(),
664 fSNNString->GetUniqueID());
666 InitializeCentBins();
668 // Print before we loop
671 // Loop over centrality bins
672 TIter next(fListOfCentralities);
673 CentralityBin* bin = 0;
674 gStyle->SetPalette(1);
675 THStack* dndetaStack = new THStack("dndeta", "dN/d#eta");
676 THStack* dndetaStackRebin = new THStack(Form("dndeta_rebin%02d", fRebin),
678 THStack* dndetaMCStack = new THStack("dndetaMC", "dN_{ch}/d#eta");
679 THStack* dndetaMCStackRebin = new THStack(Form("dndetaMC_rebin%02d", fRebin),
683 TList* truthlist = 0;
685 if (fFinalMCCorrFile.Contains(".root")) {
686 TFile* ftest = TFile::Open(fFinalMCCorrFile.Data());
688 mclist = dynamic_cast<TList*>(ftest->Get(Form("%sResults",GetName())));
689 truthlist = dynamic_cast<TList*>(ftest->Get("MCTruthResults"));
692 AliWarning("MC analysis file invalid - no final MC correction possible");
694 Int_t style = GetMarker();
695 Int_t color = GetColor();
697 AliInfo(Form("Marker style=%d, color=%d", style, color));
698 while ((bin = static_cast<CentralityBin*>(next()))) {
700 bin->End(fSums, fOutput, fNormalizationScheme, fShapeCorr, fTriggerEff,
701 fSymmetrice, fRebin, fUseROOTProj, fCorrEmpty, fCutEdges,
702 fTriggerMask, style, color, mclist, truthlist);
703 if (fCentAxis && bin->IsAllBin()) continue;
704 TH1* dndeta = bin->GetResult(0, false, "");
705 TH1* dndetaSym = bin->GetResult(0, true, "");
706 TH1* dndetaMC = bin->GetResult(0, false, "MC");
707 TH1* dndetaMCSym = bin->GetResult(0, true, "MC");
708 if (dndeta) dndetaStack->Add(dndeta);
709 if (dndetaSym) dndetaStack->Add(dndetaSym);
710 if (dndetaMC) dndetaMCStack->Add(dndetaMC);
711 if (dndetaMCSym) dndetaMCStack->Add(dndetaMCSym);
713 dndeta = bin->GetResult(fRebin, false, "");
714 dndetaSym = bin->GetResult(fRebin, true, "");
715 dndetaMC = bin->GetResult(fRebin, false, "MC");
716 dndetaMCSym = bin->GetResult(fRebin, true, "MC");
717 if (dndeta) dndetaStackRebin->Add(dndeta);
718 if (dndetaSym) dndetaStackRebin->Add(dndetaSym);
719 if (dndetaMC) dndetaMCStackRebin->Add(dndetaMC);
720 if (dndetaMCSym) dndetaMCStackRebin->Add(dndetaMCSym);
724 fOutput->Add(dndetaStack);
726 // If available output rebinned stack
727 if (!dndetaStackRebin->GetHists() ||
728 dndetaStackRebin->GetHists()->GetEntries() <= 0) {
729 AliWarning("No rebinned histograms found");
730 delete dndetaStackRebin;
731 dndetaStackRebin = 0;
733 if (dndetaStackRebin) fOutput->Add(dndetaStackRebin);
735 // If available, output track-ref stack
736 if (!dndetaMCStack->GetHists() ||
737 dndetaMCStack->GetHists()->GetEntries() <= 0) {
738 AliWarning("No MC histograms found");
739 delete dndetaMCStack;
742 if (dndetaMCStack) fOutput->Add(dndetaMCStack);
744 // If available, output rebinned track-ref stack
745 if (!dndetaMCStackRebin->GetHists() ||
746 dndetaMCStackRebin->GetHists()->GetEntries() <= 0) {
747 AliWarning("No rebinned MC histograms found");
748 delete dndetaMCStackRebin;
749 dndetaMCStackRebin = 0;
751 if (dndetaMCStackRebin) fOutput->Add(dndetaMCStackRebin);
753 // Output collision energy string
755 UShort_t sNN = fSNNString->GetUniqueID();
756 TNamed* sNNObj = new TNamed(fSNNString->GetName(),
757 AliForwardUtil::CenterOfMassEnergyString(sNN));
758 sNNObj->SetUniqueID(sNN);
759 fOutput->Add(sNNObj); // fSNNString->Clone());
762 // Output collision system string
764 UShort_t sys = fSysString->GetUniqueID();
765 TNamed* sysObj = new TNamed(fSysString->GetName(),
766 AliForwardUtil::CollisionSystemString(sys));
767 sysObj->SetUniqueID(sys);
768 fOutput->Add(sysObj); // fSysString->Clone());
771 // Output centrality axis
772 if (fCentAxis) fOutput->Add(fCentAxis);
774 // Output trigger string
775 if (fTriggerString) {
776 UShort_t mask = fTriggerString->GetUniqueID();
777 TNamed* maskObj = new TNamed(fTriggerString->GetName(),
778 AliAODForwardMult::GetTriggerString(mask));
779 maskObj->SetUniqueID(mask);
780 fOutput->Add(maskObj); // fTriggerString->Clone());
783 // Normalization string
785 UShort_t scheme = fSchemeString->GetUniqueID();
786 TNamed* schemeObj = new TNamed(fSchemeString->GetName(),
787 NormalizationSchemeString(scheme));
788 schemeObj->SetUniqueID(scheme);
789 fOutput->Add(schemeObj); // fSchemeString->Clone());
792 // Output vertex axis
793 TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax);
794 vtxAxis->SetName("vtxAxis");
795 vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
796 fOutput->Add(vtxAxis);
798 // Output trigger efficiency and shape correction
799 fOutput->Add(new TParameter<Double_t>("triggerEff", fTriggerEff));
800 if (fShapeCorr) fOutput->Add(fShapeCorr);
802 TNamed* options = new TNamed("options","");
804 str.Append(Form("Edges %scut, ", fCutEdges ? "" : "not "));
805 str.Append(Form("Empty bins %scorrected for, ", fCorrEmpty ? "" : "not "));
806 str.Append(Form("TH2::ProjectionX %sused", fUseROOTProj ? "" : "not "));
807 options->SetTitle(str);
808 fOutput->Add(options);
810 PostData(2, fOutput);
812 //________________________________________________________________________
814 AliBasedNdetaTask::LoadNormalizationData(UShort_t sys, UShort_t energy)
816 // Load the normalisation data for dN/deta for pp INEL, INEL>0, and NSD
817 DGUARD(fDebug,1,"Load the normalization data for sys=%d, energy=%d",
821 if(energy == 7000) snn.Form("7000");
822 if(energy == 2750) snn.Form("2750");
824 if(fShapeCorr && (fTriggerEff != 1)) {
825 AliInfo("Objects already set for normalization - no action taken");
829 TFile* fin = TFile::Open(Form("$ALICE_ROOT/PWGLF/FORWARD/corrections/"
830 "Normalization/normalizationHists_%s_%s.root",
831 type.Data(),snn.Data()));
833 AliWarning(Form("no file for normalization of %d/%d", sys, energy));
838 if ((fNormalizationScheme & kShape) && !fShapeCorr) {
839 TString trigName("All");
840 if (fTriggerMask == AliAODForwardMult::kInel ||
841 fTriggerMask == AliAODForwardMult::kNClusterGt0)
843 else if (fTriggerMask == AliAODForwardMult::kNSD)
845 else if (fTriggerMask == AliAODForwardMult::kInelGt0)
846 trigName = "InelGt0";
848 AliWarning(Form("Normalization for trigger %s not known, using all",
849 AliAODForwardMult::GetTriggerString(fTriggerMask)));
852 TString shapeCorName(Form("h%sNormalization", trigName.Data()));
853 TH2F* shapeCor = dynamic_cast<TH2F*>(fin->Get(shapeCorName));
854 if (shapeCor) SetShapeCorrection(shapeCor);
856 AliWarning(Form("No shape correction found for %s", trigName.Data()));
860 // Trigger efficiency
861 TString effName(Form("%sTriggerEff",
862 fTriggerMask == AliAODForwardMult::kInel ? "inel" :
863 fTriggerMask == AliAODForwardMult::kNSD ? "nsd" :
864 fTriggerMask == AliAODForwardMult::kInelGt0 ?
866 TParameter<float>* eff = 0;
867 if (fNormalizationScheme & kTriggerEfficiency)
868 eff = static_cast<TParameter<float>*>(fin->Get(effName));
869 Double_t trigEff = eff ? eff->GetVal() : 1;
870 if (fTriggerEff != 1) SetTriggerEff(trigEff);
871 if (fTriggerEff < 0) fTriggerEff = 1;
874 // Rescale the shape correction by the trigger efficiency
876 AliWarning(Form("Rescaling shape correction by trigger efficiency: "
877 "1/E_X=1/%f", trigEff));
878 fShapeCorr->Scale(1. / trigEff);
882 if (fShapeCorr && fTriggerEff) AliInfo("Loaded objects for normalization.");
886 //________________________________________________________________________
888 AliBasedNdetaTask::Print(Option_t*) const
893 std::cout << this->ClassName() << ": " << this->GetName() << "\n"
895 << " Trigger: " << (fTriggerString ?
896 fTriggerString->GetTitle() :
898 << " Vertex range: [" << fVtxMin << ":" << fVtxMax << "]\n"
899 << " Rebin factor: " << fRebin << "\n"
900 << " Cut edges: " << fCutEdges << "\n"
901 << " Symmertrice: " << fSymmetrice << "\n"
902 << " Use TH2::ProjectionX: " << fUseROOTProj << "\n"
903 << " Correct for empty: " << fCorrEmpty << "\n"
904 << " Normalization scheme: " << (fSchemeString ?
905 fSchemeString->GetTitle() :
907 << " Trigger efficiency: " << fTriggerEff << "\n"
908 << " Shape correction: " << (fShapeCorr ?
909 fShapeCorr->GetName() :
911 << " sqrt(s_NN): " << (fSNNString ?
912 fSNNString->GetTitle() :
914 << " Collision system: " << (fSysString ?
915 fSysString->GetTitle() :
917 << " Centrality bins: " << (fCentAxis ? "" : "none");
919 Int_t nBins = fCentAxis->GetNbins();
920 const Double_t* bins = fCentAxis->GetXbins()->GetArray();
921 for (Int_t i = 0; i <= nBins; i++)
922 std::cout << (i==0 ? "" : "-") << bins[i];
924 std::cout << std::noboolalpha << std::endl;
928 //________________________________________________________________________
930 AliBasedNdetaTask::Rebin(const TH1D* h, Int_t rebin, Bool_t cutEdges)
933 // Make a copy of the input histogram and rebin that histogram
936 // h Histogram to rebin
939 // New (rebinned) histogram
941 if (rebin <= 1) return 0;
943 Int_t nBins = h->GetNbinsX();
944 if(nBins % rebin != 0) {
945 AliWarningGeneral("AliBasedNdetaTask",
946 Form("Rebin factor %d is not a devisor of current number "
947 "of bins %d in the histogram %s",
948 rebin, nBins, h->GetName()));
953 TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d",
954 h->GetName(), rebin)));
956 tmp->SetDirectory(0);
958 // The new number of bins
959 Int_t nBinsNew = nBins / rebin;
960 for(Int_t i = 1;i<= nBinsNew; i++) {
961 Double_t content = 0;
965 for(Int_t j = 1; j<=rebin;j++) {
966 Int_t bin = (i-1)*rebin + j;
967 Double_t c = h->GetBinContent(bin);
968 if (c <= 0) continue;
971 if (h->GetBinContent(bin+1)<=0 ||
972 h->GetBinContent(bin-1)<=0) {
974 AliWarningGeneral("AliBasedNdetaTask",
975 Form("removing bin %d=%f of %s (%d=%f,%d=%f)",
976 bin, c, h->GetName(),
977 bin+1, h->GetBinContent(bin+1),
978 bin-1, h->GetBinContent(bin-1)));
983 Double_t e = h->GetBinError(bin);
984 Double_t w = 1 / (e*e); // 1/c/c
991 if(content > 0 && nbins > 0) {
992 tmp->SetBinContent(i, wsum / sumw);
993 tmp->SetBinError(i,1./TMath::Sqrt(sumw));
1000 //__________________________________________________________________
1002 AliBasedNdetaTask::Symmetrice(const TH1* h)
1005 // Make an extension of @a h to make it symmetric about 0
1008 // h Histogram to symmertrice
1011 // Symmetric extension of @a h
1013 Int_t nBins = h->GetNbinsX();
1014 TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
1015 s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
1017 s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
1018 s->SetMarkerColor(h->GetMarkerColor());
1019 s->SetMarkerSize(h->GetMarkerSize());
1020 s->SetMarkerStyle(FlipHollowStyle(h->GetMarkerStyle()));
1021 s->SetFillColor(h->GetFillColor());
1022 s->SetFillStyle(h->GetFillStyle());
1025 // Find the first and last bin with data
1026 Int_t first = nBins+1;
1028 for (Int_t i = 1; i <= nBins; i++) {
1029 if (h->GetBinContent(i) <= 0) continue;
1030 first = TMath::Min(first, i);
1031 last = TMath::Max(last, i);
1034 Double_t xfirst = h->GetBinCenter(first-1);
1035 Int_t f1 = h->GetXaxis()->FindBin(-xfirst);
1036 Int_t l2 = s->GetXaxis()->FindBin(xfirst);
1037 for (Int_t i = f1, j=l2; i <= last; i++,j--) {
1038 s->SetBinContent(j, h->GetBinContent(i));
1039 s->SetBinError(j, h->GetBinError(i));
1041 // Fill in overlap bin
1042 s->SetBinContent(l2+1, h->GetBinContent(first));
1043 s->SetBinError(l2+1, h->GetBinError(first));
1047 //__________________________________________________________________
1049 AliBasedNdetaTask::GetMarkerStyle(UShort_t bits)
1051 Int_t base = bits & (0xFE);
1052 Bool_t hollow = bits & kHollow;
1054 case kCircle: return (hollow ? 24 : 20);
1055 case kSquare: return (hollow ? 25 : 21);
1056 case kUpTriangle: return (hollow ? 26 : 22);
1057 case kDownTriangle: return (hollow ? 32 : 23);
1058 case kDiamond: return (hollow ? 27 : 33);
1059 case kCross: return (hollow ? 28 : 34);
1060 case kStar: return (hollow ? 30 : 29);
1064 //__________________________________________________________________
1066 AliBasedNdetaTask::GetMarkerBits(Int_t style)
1070 case 24: case 25: case 26: case 27: case 28: case 30: case 32:
1071 bits |= kHollow; break;
1074 case 20: case 24: bits |= kCircle; break;
1075 case 21: case 25: bits |= kSquare; break;
1076 case 22: case 26: bits |= kUpTriangle; break;
1077 case 23: case 32: bits |= kDownTriangle; break;
1078 case 27: case 33: bits |= kDiamond; break;
1079 case 28: case 34: bits |= kCross; break;
1080 case 29: case 30: bits |= kStar; break;
1084 //__________________________________________________________________
1086 AliBasedNdetaTask::FlipHollowStyle(Int_t style)
1088 UShort_t bits = GetMarkerBits(style);
1089 Int_t ret = GetMarkerStyle(bits ^ kHollow);
1093 //====================================================================
1095 AliBasedNdetaTask::Sum::Init(TList* list, const TH2D* data, Int_t col)
1097 DGUARD(fDebug,1,"Initializing sums with %s", data->GetName());
1098 TString n(GetHistName(0));
1099 TString n0(GetHistName(1));
1100 const char* postfix = GetTitle();
1102 fSum = static_cast<TH2D*>(data->Clone(n));
1103 if (postfix) fSum->SetTitle(Form("%s (%s)", data->GetTitle(), postfix));
1104 fSum->SetDirectory(0);
1105 fSum->SetMarkerColor(col);
1106 fSum->SetMarkerStyle(GetMarkerStyle(kCircle|kSolid));
1110 fSum0 = static_cast<TH2D*>(data->Clone(n0));
1112 fSum0->SetTitle(Form("%s 0-bin (%s)", data->GetTitle(), postfix));
1114 fSum0->SetTitle(Form("%s 0-bin", data->GetTitle()));
1115 fSum0->SetDirectory(0);
1116 fSum0->SetMarkerColor(col);
1117 fSum0->SetMarkerStyle(GetMarkerStyle(kCross|kHollow));
1121 fEvents = new TH1I(GetHistName(2), "Event types", 2, -.5, 1.5);
1122 fEvents->SetDirectory(0);
1123 fEvents->GetXaxis()->SetBinLabel(1, "Non-zero");
1124 fEvents->GetXaxis()->SetBinLabel(2, "Zero");
1128 //____________________________________________________________________
1130 AliBasedNdetaTask::Sum::GetHistName(const char* name, Int_t what, const char* post)
1133 if (what == 1) n.Append("0");
1134 else if (what == 2) n.Append("Events");
1135 if (post && post[0] != '\0') n.Append(post);
1139 //____________________________________________________________________
1141 AliBasedNdetaTask::Sum::GetHistName(Int_t what) const
1143 return GetHistName(GetName(), what, GetTitle());
1146 //____________________________________________________________________
1148 AliBasedNdetaTask::Sum::Add(const TH2D* data, Bool_t isZero)
1150 DGUARD(fDebug,2,"Adding %s to sums", data->GetName());
1151 if (isZero) fSum0->Add(data);
1152 else fSum->Add(data);
1153 fEvents->Fill(isZero ? 1 : 0);
1156 //____________________________________________________________________
1158 AliBasedNdetaTask::Sum::CalcSum(TList* output,
1164 Bool_t corrEmpty) const
1166 DGUARD(fDebug,2,"Calculating final summed histogram %s", fSum->GetName());
1168 TH2D* ret = static_cast<TH2D*>(fSum->Clone(fSum->GetName()));
1169 ret->SetDirectory(0);
1171 Int_t n = Int_t(fEvents->GetBinContent(1));
1172 Int_t n0 = Int_t(fEvents->GetBinContent(2));
1174 AliInfoF("Adding histograms %s and %s with weights %f and %f resp.",
1175 fSum0->GetName(), fSum->GetName(), 1./epsilon, 1./epsilon0);
1176 DMSG(fDebug,2, "Adding histograms %s and %s with weights %f and %f resp.",
1177 fSum0->GetName(), fSum->GetName(), 1./epsilon, 1./epsilon0);
1178 // Generate merged histogram
1179 ret->Add(fSum0, fSum, 1. / epsilon0, 1. / epsilon);
1180 ntotal = n / epsilon + n0 / epsilon0;
1182 TList* out = new TList;
1184 const char* postfix = GetTitle();
1185 if (!postfix) postfix = "";
1186 out->SetName(Form("partial%s", postfix));
1189 // Now make copies, normalize them, and store in output list
1190 TH2D* sumCopy = static_cast<TH2D*>(fSum->Clone("sum"));
1191 TH2D* sum0Copy = static_cast<TH2D*>(fSum0->Clone("sum0"));
1192 TH2D* retCopy = static_cast<TH2D*>(ret->Clone("sumAll"));
1193 sumCopy->SetMarkerStyle(FlipHollowStyle(marker));
1194 sumCopy->SetDirectory(0);
1195 sum0Copy->SetMarkerStyle(GetMarkerStyle(GetMarkerBits(marker)+4));
1196 sum0Copy->SetDirectory(0);
1197 retCopy->SetMarkerStyle(marker);
1198 retCopy->SetDirectory(0);
1200 TH1D* norm = ProjectX(fSum, "norm", 0, 0, rootProj, corrEmpty, false);
1201 TH1D* norm0 = ProjectX(fSum0, "norm0", 0, 0, rootProj, corrEmpty, false);
1202 TH1D* normAll = ProjectX(ret, "normAll", 0, 0, rootProj, corrEmpty, false);
1203 norm->SetDirectory(0);
1204 norm0->SetDirectory(0);
1205 normAll->SetDirectory(0);
1207 ScaleToCoverage(sumCopy, norm);
1208 ScaleToCoverage(sum0Copy, norm0);
1209 ScaleToCoverage(retCopy, normAll);
1211 Int_t nY = fSum->GetNbinsY();
1212 TH1D* sumCopyPx = ProjectX(sumCopy, "average", 1, nY,rootProj,corrEmpty);
1213 TH1D* sum0CopyPx = ProjectX(sum0Copy, "average0", 1, nY,rootProj,corrEmpty);
1214 TH1D* retCopyPx = ProjectX(retCopy, "averageAll", 1, nY,rootProj,corrEmpty);
1215 sumCopyPx->SetDirectory(0);
1216 sum0CopyPx->SetDirectory(0);
1217 retCopyPx->SetDirectory(0);
1219 // Scale our 1D histograms
1220 sumCopyPx->Scale(1., "width");
1221 sum0CopyPx->Scale(1., "width");
1222 retCopyPx->Scale(1., "width");
1224 AliInfo(Form("Maximum %f,%f changed to %f", sumCopyPx->GetMaximum(),
1225 sum0CopyPx->GetMaximum(), retCopyPx->GetMaximum()));
1227 // Scale the normalization - they should be 1 at the maximum
1228 norm->Scale(n > 0 ? 1. / n : 1);
1229 norm0->Scale(n0 > 0 ? 1. / n0 : 1);
1230 normAll->Scale(ntotal > 0 ? 1. / ntotal : 1);
1235 out->Add(sumCopyPx);
1236 out->Add(sum0CopyPx);
1237 out->Add(retCopyPx);
1242 AliInfoF("Returning (1/%f * %s + 1/%f * %s), "
1243 "1/%f * %d + 1/%f * %d = %d",
1244 epsilon0, fSum0->GetName(), epsilon, fSum->GetName(),
1245 epsilon0, n0, epsilon, n, int(ntotal));
1247 for (Int_t i = 1; i <= ret->GetNbinsX(); i++) {
1248 Double_t nc = sum->GetBinContent(i, 0);
1249 Double_t nc0 = sum0->GetBinContent(i, 0);
1250 ret->SetBinContent(i, 0, nc + nc0); // Just count events
1257 //====================================================================
1258 AliBasedNdetaTask::CentralityBin::CentralityBin()
1267 fDoFinalMCCorrection(false),
1273 DGUARD(0,0,"Centrality bin default construction");
1275 //____________________________________________________________________
1276 AliBasedNdetaTask::CentralityBin::CentralityBin(const char* name,
1277 Short_t low, Short_t high)
1286 fDoFinalMCCorrection(false),
1293 // name Name used for histograms (e.g., Forward)
1294 // low Lower centrality cut in percent
1295 // high Upper centrality cut in percent
1297 DGUARD(0,0,"Named Centrality bin construction: %s [%d,%d]", name);
1298 if (low <= 0 && high <= 0) {
1301 SetTitle("All centralities");
1306 SetTitle(Form("Centrality bin from %3d%% to %3d%%", low, high));
1309 //____________________________________________________________________
1310 AliBasedNdetaTask::CentralityBin::CentralityBin(const CentralityBin& o)
1316 fTriggers(o.fTriggers),
1319 fDoFinalMCCorrection(o.fDoFinalMCCorrection),
1326 // other Object to copy from
1328 DGUARD(0,0,"Copy Centrality bin construction");
1330 //____________________________________________________________________
1331 AliBasedNdetaTask::CentralityBin::~CentralityBin()
1336 DGUARD(fDebug,3,"Centrality bin desctruction");
1337 if (fSums) fSums->Delete();
1338 if (fOutput) fOutput->Delete();
1341 //____________________________________________________________________
1342 AliBasedNdetaTask::CentralityBin&
1343 AliBasedNdetaTask::CentralityBin::operator=(const CentralityBin& o)
1346 // Assignment operator
1349 // other Object to assign from
1352 // Reference to this
1354 DGUARD(fDebug,3,"Centrality bin assignment");
1355 if (&o == this) return *this;
1356 SetName(o.GetName());
1357 SetTitle(o.GetTitle());
1359 fOutput = o.fOutput;
1362 fTriggers = o.fTriggers;
1365 fDoFinalMCCorrection = o.fDoFinalMCCorrection;
1369 //____________________________________________________________________
1371 AliBasedNdetaTask::CentralityBin::GetColor(Int_t fallback) const
1373 if (IsAllBin()) return fallback;
1374 Float_t fc = (fLow+double(fHigh-fLow)/2) / 100;
1375 Int_t nCol = gStyle->GetNumberOfColors();
1376 Int_t icol = TMath::Min(nCol-1,int(fc * nCol + .5));
1377 Int_t col = gStyle->GetColorPalette(icol);
1380 //____________________________________________________________________
1382 AliBasedNdetaTask::CentralityBin::GetListName() const
1385 // Get the list name
1390 if (IsAllBin()) return "all";
1391 return Form("cent%03d_%03d", fLow, fHigh);
1393 //____________________________________________________________________
1395 AliBasedNdetaTask::CentralityBin::CreateOutputObjects(TList* dir)
1398 // Create output objects
1403 DGUARD(fDebug,1,"Create centrality bin output objects");
1405 fSums->SetName(GetListName());
1409 fTriggers = AliAODForwardMult::MakeTriggerHistogram("triggers");
1410 fTriggers->SetDirectory(0);
1411 fSums->Add(fTriggers);
1413 //____________________________________________________________________
1415 AliBasedNdetaTask::CentralityBin::SetDebugLevel(Int_t lvl)
1418 if (fSum) fSum->fDebug = lvl;
1419 if (fSumMC) fSumMC->fDebug = lvl;
1422 //____________________________________________________________________
1424 AliBasedNdetaTask::CentralityBin::ReadSum(TList* list, bool mc)
1426 const char* post = (mc ? "MC" : "");
1427 TString sn = Sum::GetHistName(GetName(),0,post);
1428 TString sn0 = Sum::GetHistName(GetName(),1,post);
1429 TString ev = Sum::GetHistName(GetName(),2,post);
1430 TH2D* sum = static_cast<TH2D*>(list->FindObject(sn));
1431 TH2D* sum0 = static_cast<TH2D*>(list->FindObject(sn0));
1432 TH1I* events = static_cast<TH1I*>(list->FindObject(ev));
1433 if (!sum || !sum0 || !events) {
1434 AliWarningF("Failed to find one or more histograms: "
1435 "%s (%p) %s (%p) %s (%p)",
1441 Sum* ret = new Sum(GetName(), post);
1444 ret->fEvents = events;
1445 ret->fDebug = fDebug;
1446 if (mc) fSumMC = ret;
1452 //____________________________________________________________________
1454 AliBasedNdetaTask::CentralityBin::CreateSums(const TH2D* data, const TH2D* mc)
1457 // Create sum histogram
1460 // data Data histogram to clone
1461 // mc (optional) MC histogram to clone
1463 DGUARD(fDebug,1,"Create centrality bin sums from %s",
1464 data ? data->GetName() : "(null)");
1466 fSum = new Sum(GetName(),"");
1467 fSum->Init(fSums, data, GetColor());
1468 fSum->fDebug = fDebug;
1471 // If no MC data is given, then do not create MC sum histogram
1474 fSumMC = new Sum(GetName(), "MC");
1475 fSumMC->Init(fSums, mc, GetColor());
1476 fSumMC->fDebug = fDebug;
1479 //____________________________________________________________________
1481 AliBasedNdetaTask::CentralityBin::CheckEvent(const AliAODForwardMult* forward,
1483 Double_t vzMin, Double_t vzMax)
1486 // Check the trigger, vertex, and centrality
1492 // true if the event is to be used
1494 if (!forward) return false;
1496 DGUARD(fDebug,2,"Check the event");
1497 // We do not check for centrality here - it's already done
1498 return forward->CheckEvent(triggerMask, vzMin, vzMax, 0, 0, fTriggers);
1502 //____________________________________________________________________
1504 AliBasedNdetaTask::CentralityBin::ProcessEvent(const AliAODForwardMult* forward,
1505 Int_t triggerMask, Bool_t isZero,
1506 Double_t vzMin, Double_t vzMax,
1507 const TH2D* data, const TH2D* mc)
1513 // forward Forward data (for trigger, vertex, & centrality)
1514 // triggerMask Trigger mask
1515 // vzMin Minimum IP z coordinate
1516 // vzMax Maximum IP z coordinate
1517 // data Data histogram
1520 DGUARD(fDebug,1,"Process one event for %s a given centrality bin",
1521 data ? data->GetName() : "(null)");
1522 if (!CheckEvent(forward, triggerMask, vzMin, vzMax)) return;
1524 if (!fSum) CreateSums(data, mc);
1526 fSum->Add(data, isZero);
1527 if (mc) fSumMC->Add(mc, isZero);
1530 //________________________________________________________________________
1532 AliBasedNdetaTask::CentralityBin::Normalization(const TH1I& t,
1535 Double_t& ntotal) const
1538 // Calculate normalization
1541 // t Trigger histogram
1542 // scheme Normaliztion scheme
1544 // ntotal On return, contains the number of events.
1546 DGUARD(fDebug,1,"Normalize centrality bin with %s", t.GetName());
1547 Double_t nAll = t.GetBinContent(AliAODForwardMult::kBinAll);
1548 Double_t nB = t.GetBinContent(AliAODForwardMult::kBinB);
1549 Double_t nA = t.GetBinContent(AliAODForwardMult::kBinA);
1550 Double_t nC = t.GetBinContent(AliAODForwardMult::kBinC);
1551 Double_t nE = t.GetBinContent(AliAODForwardMult::kBinE);
1552 Double_t nOffline = t.GetBinContent(AliAODForwardMult::kBinOffline);
1553 Double_t nTriggered = t.GetBinContent(AliAODForwardMult::kWithTrigger);
1554 Double_t nWithVertex = t.GetBinContent(AliAODForwardMult::kWithVertex);
1555 Double_t nAccepted = ntotal; // t.GetBinContent(AliAODForwardMult::kAccepted);
1558 if (nTriggered <= 0.1) {
1559 AliError("Number of triggered events <= 0");
1562 if (nWithVertex <= 0.1) {
1563 AliError("Number of events with vertex <= 0");
1567 Double_t vtxEff = nWithVertex / nTriggered;
1568 Double_t scaler = 1;
1569 Double_t beta = nA + nC - 2*nE;
1571 if (scheme & kEventLevel && !(scheme & kZeroBin)) {
1572 ntotal = nAccepted / vtxEff;
1574 AliInfo(Form("Calculating event normalisation as\n"
1575 " N = N_A * N_T / N_V = %d * %d / %d = %f (%f)",
1576 Int_t(nAccepted), Int_t(nTriggered), Int_t(nWithVertex),
1579 if (scheme & kBackground) {
1581 // s = --------- = ------------- = ------------
1582 // 1 - beta 1 - beta E_V 1 - beta N_V
1583 // --- ---- -------- ---- ---
1584 // E_V N_V N_V N_V N_T
1592 ntotal -= nAccepted * beta / nWithVertex;
1593 // This one is direct and correct.
1594 // scaler = 1. / (1. / vtxEff - beta / nWithVertex);
1595 // A simpler expresion
1596 scaler /= (1 - beta / nTriggered); // 0.831631 -> 0.780689
1597 AliInfo(Form("Calculating event normalisation as\n"
1598 " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n"
1599 " N = N - N_A * beta / N_V = %f - %d * %d / %d = %f (%f)",
1600 Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta),
1601 nAccepted / vtxEff, Int_t(nAccepted), Int_t(beta),
1602 Int_t(nWithVertex), ntotal, scaler));
1605 if (scheme & kZeroBin) {
1608 // N = N_A + 1/E_X * N_A / N_V (N_T - N_V - beta)
1609 // = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V))
1611 // s = N_A/N = 1 / (1 + 1/E_X (N_T/N_V - 1 - beta / N_V))
1612 // = N_V / (N_V + 1/E_X (N_T - N_V - beta))
1614 if (!(scheme & kBackground)) beta = 0;
1615 ntotal = nAccepted * (1 + 1/trigEff * (nTriggered / nWithVertex - 1
1616 - beta / nWithVertex));
1617 scaler = nWithVertex / (nWithVertex +
1618 1/trigEff * (nTriggered-nWithVertex-beta));
1619 AliInfo(Form("Calculating event normalisation as\n"
1620 " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n"
1621 " N = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V)) = "
1622 "%d (1 + 1 / %f (%d / %d - 1 - %d / %d)) = %f (%f)",
1623 Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta),
1624 Int_t(nAccepted), trigEff, Int_t(nTriggered),
1625 Int_t(nWithVertex), Int_t(beta), Int_t(nWithVertex),
1628 if (scheme & kTriggerEfficiency && !(scheme & kZeroBin)) {
1631 AliInfo(Form("Correcting for trigger efficiency:\n"
1632 " N = 1 / E_X * N = 1 / %f * %d = %f (%f)",
1633 trigEff, Int_t(ntotal), ntotal / trigEff, scaler));
1637 " Total of %9d events for %s\n"
1638 " of these %9d have an offline trigger\n"
1639 " of these N_T = %9d has the selected trigger\n"
1640 " of these N_V = %9d has a vertex\n"
1641 " of these N_A = %9d were in the selected range\n"
1642 " Triggers by hardware type:\n"
1644 " N_ac = %9d (%d+%d)\n"
1646 " Vertex efficiency: %f\n"
1647 " Trigger efficiency: %f\n"
1648 " Total number of events: N = %f\n"
1649 " Scaler (N_A/N): %f",
1650 Int_t(nAll), GetTitle(), Int_t(nOffline),
1651 Int_t(nTriggered), Int_t(nWithVertex), Int_t(nAccepted),
1652 Int_t(nB), Int_t(nA+nC), Int_t(nA), Int_t(nC), Int_t(nE),
1653 vtxEff, trigEff, ntotal, scaler));
1657 //________________________________________________________________________
1659 AliBasedNdetaTask::CentralityBin::GetResultName(Int_t rebin,
1661 const char* postfix) const
1664 n = Form("dndeta%s%s",GetName(), postfix);
1665 if (rebin > 1) n.Append(Form("_rebin%02d", rebin));
1666 if (sym) n.Append("_mirror");
1669 //________________________________________________________________________
1671 AliBasedNdetaTask::CentralityBin::GetResult(Int_t rebin,
1673 const char* postfix) const
1676 AliWarning(Form("No output list defined in %s [%3d,%3d]", GetName(),
1680 TString n = GetResultName(rebin, sym, postfix);
1681 TObject* o = fOutput->FindObject(n.Data());
1683 // AliWarning(Form("Object %s not found in output list", n.Data()));
1686 return static_cast<TH1*>(o);
1689 //________________________________________________________________________
1691 AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum,
1692 const char* postfix,
1695 const TH2F* shapeCorr,
1706 // Generate the dN/deta result from input
1709 // sum Sum of 2D hists
1710 // postfix Post fix on names
1711 // rootProj Whether to use ROOT TH2::ProjectionX
1712 // corrEmpty Correct for empty bins
1713 // shapeCorr Shape correction to use
1714 // scaler Event-level normalization scaler
1715 // symmetrice Whether to make symmetric extensions
1716 // rebin Whether to rebin
1717 // cutEdges Whether to cut edges when rebinning
1719 DGUARD(fDebug,1,"Make centrality bin result from %s", sum->GetName());
1720 TH2D* copy = static_cast<TH2D*>(sum->Clone(Form("d2Ndetadphi%s%s",
1721 GetName(), postfix)));
1722 TH1D* accNorm = ProjectX(sum, Form("norm%s%s",GetName(), postfix), 0, 0,
1723 rootProj, corrEmpty, false);
1724 accNorm->SetDirectory(0);
1726 // ---- Scale by shape correction ----------------------------------
1727 if (shapeCorr) copy->Divide(shapeCorr);
1728 else AliInfo("No shape correction specified, or disabled");
1730 // --- Normalize to the coverage -----------------------------------
1731 ScaleToCoverage(copy, accNorm);
1733 // --- Event-level normalization -----------------------------------
1734 copy->Scale(scaler);
1736 // --- Project on X axis -------------------------------------------
1737 TH1D* dndeta = ProjectX(copy, Form("dndeta%s%s",GetName(), postfix),
1738 1, sum->GetNbinsY(), rootProj, corrEmpty);
1739 dndeta->SetDirectory(0);
1740 // Event-level normalization
1741 dndeta->Scale(1., "width");
1742 copy->Scale(1., "width");
1744 TH1D* dndetaMCCorrection = 0;
1745 TList* centlist = 0;
1746 TH1D* dndetaMCtruth = 0;
1747 TList* truthcentlist = 0;
1749 // Possible final correction to <MC analysis> / <MC truth>
1751 centlist = static_cast<TList*> (mclist->FindObject(GetListName()));
1753 dndetaMCCorrection = static_cast<TH1D*>(centlist->FindObject(Form("dndeta%s%s",GetName(), postfix)));
1755 truthcentlist = static_cast<TList*> (truthlist->FindObject(GetListName()));
1757 dndetaMCtruth = static_cast<TH1D*> (truthcentlist->FindObject("dndetaTruth"));
1758 //std::cout<<dndetaMCCorrection<<" "<<dndetaMCtruth<<std::endl;
1759 if(dndetaMCCorrection && dndetaMCtruth) {
1760 AliInfo("Correcting with final MC correction");
1761 dndetaMCCorrection->Divide(dndetaMCtruth);
1762 dndeta->Divide(dndetaMCCorrection);
1764 //std::cout<<"histo "<<Form("dndeta%s%s",GetName(), postfix)<<" "<<GetListName()<<" "<<dndetaMCCorrection<<std::endl;
1765 //std::cout<<"truth "<<GetListName()<<" "<<dndetaMCtruth<<std::endl;
1768 else AliInfo("No final MC correction applied");
1770 // --- Set some histogram attributes -------------------------------
1772 Int_t rColor = GetColor(color);
1773 if (postfix && postfix[0] != '\0') post = Form(" (%s)", postfix);
1774 SetHistogramAttributes(dndeta, rColor, marker,
1775 Form("ALICE %s%s", GetName(), post.Data()));
1776 SetHistogramAttributes(accNorm, rColor, marker,
1777 Form("ALICE %s normalisation%s",
1778 GetName(), post.Data()));
1780 // --- Make symmetric extensions and rebinnings --------------------
1781 if (symmetrice) fOutput->Add(Symmetrice(dndeta));
1782 fOutput->Add(dndeta);
1783 fOutput->Add(accNorm);
1785 fOutput->Add(Rebin(dndeta, rebin, cutEdges));
1786 if (symmetrice) fOutput->Add(Symmetrice(Rebin(dndeta, rebin, cutEdges)));
1789 //________________________________________________________________________
1791 AliBasedNdetaTask::CentralityBin::End(TList* sums,
1794 const TH2F* shapeCorr,
1808 // End of processing
1811 // sums List of sums
1812 // results Output list of results
1813 // shapeCorr Shape correction or nil
1814 // trigEff Trigger efficiency
1815 // symmetrice Whether to symmetrice the results
1816 // rebin Whether to rebin the results
1817 // corrEmpty Whether to correct for empty bins
1818 // cutEdges Whether to cut edges when rebinning
1819 // triggerMask Trigger mask
1821 DGUARD(fDebug,1,"End centrality bin procesing");
1823 fSums = dynamic_cast<TList*>(sums->FindObject(GetListName()));
1825 AliError("Could not retrieve TList fSums");
1829 fOutput = new TList;
1830 fOutput->SetName(GetListName());
1831 fOutput->SetOwner();
1832 results->Add(fOutput);
1835 if (!ReadSum(fSums, false)) {
1836 AliInfo("This task did not produce any output");
1840 if (!fSumMC) ReadSum(fSums, true);
1842 fTriggers = static_cast<TH1I*>(fSums->FindObject("triggers"));
1845 AliError("Couldn't find histogram 'triggers' in list");
1849 // --- Get normalization scaler ------------------------------------
1850 Double_t epsilonT = trigEff;
1851 Double_t epsilonT0 = trigEff;
1853 if (triggerMask == AliAODForwardMult::kNSD) {
1854 // This is a local change
1856 AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT));
1858 else if (triggerMask == AliAODForwardMult::kInel) {
1859 // This is a local change
1861 AliWarning(Form("Using hard-coded Inel trigger efficiency of %f",epsilonT));
1863 if (scheme & kZeroBin) {
1864 if (triggerMask==AliAODForwardMult::kInel)
1865 epsilonT0 = 0.785021; // 0.100240;
1866 else if (triggerMask==AliAODForwardMult::kInelGt0)
1868 else if (triggerMask==AliAODForwardMult::kNSD)
1869 epsilonT0 = .706587;
1871 AliWarning(Form("Using hard-coded NCluster>0 trigger efficiency of %f",
1875 // Get our histograms
1877 TH2D* sum = fSum->CalcSum(fOutput, nSum, epsilonT0, 1,
1878 marker, rootProj, corrEmpty);
1879 Double_t nSumMC = 0;
1881 if (fSumMC) sumMC = fSumMC->CalcSum(fOutput, nSumMC,
1882 epsilonT0, 1, marker,
1883 rootProj, corrEmpty);
1885 AliError("Failed to get sum from summer - bailing out");
1889 Double_t ntotal = nSum;
1890 Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal);
1892 AliError("Failed to calculate normalization - bailing out");
1895 fOutput->Add(fTriggers->Clone());
1897 // --- Make result and store ---------------------------------------
1898 MakeResult(sum, "", rootProj, corrEmpty, (scheme & kShape) ? shapeCorr : 0,
1899 scaler, symmetrice, rebin, cutEdges, marker, color,
1902 // --- Process result from TrackRefs -------------------------------
1904 MakeResult(sumMC, "MC", rootProj, corrEmpty,
1905 (scheme & kShape) ? shapeCorr : 0,
1906 scaler, symmetrice, rebin, cutEdges,
1907 GetMarkerStyle(GetMarkerBits(marker)+4), color,
1911 // if (!IsAllBin()) return;