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);
202 //________________________________________________________________________
204 AliBasedNdetaTask::SetNormalizationScheme(const char* what)
207 // Set normalisation scheme
209 DGUARD(fDebug,3,"Set the normalization scheme: %s", what);
214 TIter next(twhat.Tokenize(" ,|"));
215 while ((opt = static_cast<TObjString*>(next()))) {
216 TString s(opt->GetString());
217 if (s.IsNull()) continue;
220 case '-': add = false; // Fall through
221 case '+': s.Remove(0,1); // Remove character
224 if (s.CompareTo("EVENT") == 0) bit = kEventLevel;
225 else if (s.CompareTo("SHAPE") == 0) bit = kShape;
226 else if (s.CompareTo("BACKGROUND")== 0) bit = kBackground;
227 else if (s.CompareTo("TRIGGER") == 0) bit = kTriggerEfficiency;
228 else if (s.CompareTo("FULL") == 0) bit = kFull;
229 else if (s.CompareTo("NONE") == 0) bit = kNone;
230 else if (s.CompareTo("ZEROBIN") == 0) bit = kZeroBin;
232 Warning("SetNormalizationScheme", "Unknown option %s", s.Data());
233 if (add) scheme |= bit;
236 SetNormalizationScheme(scheme);
238 //________________________________________________________________________
240 AliBasedNdetaTask::SetNormalizationScheme(UShort_t scheme)
242 DGUARD(fDebug,3,"Set the normalization scheme: 0x%x", scheme);
243 fNormalizationScheme = scheme;
245 if (scheme == kFull) tit = "FULL";
247 if (scheme & kEventLevel) tit.Append("EVENT ");
248 if (scheme & kShape) tit.Append("SHAPE ");
249 if (scheme & kBackground) tit.Append("BACKGROUND ");
250 if (scheme & kTriggerEfficiency) tit.Append("TRIGGER ");
251 if (scheme & kZeroBin) tit.Append("ZEROBIN ");
253 tit = tit.Strip(TString::kBoth);
254 if (!fSchemeString) fSchemeString = new TNamed("scheme", "");
255 fSchemeString->SetTitle(tit);
256 fSchemeString->SetUniqueID(fNormalizationScheme);
258 //________________________________________________________________________
260 AliBasedNdetaTask::SetTriggerMask(const char* mask)
263 // Set the trigger maskl
268 DGUARD(fDebug,3,"Set the trigger mask: %s", mask);
269 SetTriggerMask(AliAODForwardMult::MakeTriggerMask(mask));
271 //________________________________________________________________________
273 AliBasedNdetaTask::SetTriggerMask(UShort_t mask)
275 DGUARD(fDebug,3,"Set the trigger mask: 0x%0x", mask);
277 TString tit(AliAODForwardMult::GetTriggerString(mask));
278 tit = tit.Strip(TString::kBoth);
279 if (!fTriggerString) fTriggerString = new TNamed("trigger", "");
280 fTriggerString->SetTitle(tit);
281 fTriggerString->SetUniqueID(fTriggerMask);
284 //________________________________________________________________________
286 AliBasedNdetaTask::SetShapeCorrection(const TH2F* c)
289 // Set the shape correction (a.k.a., track correction) for selected
295 DGUARD(fDebug,3,"Set the shape correction: %p", c);
297 fShapeCorr = static_cast<TH2F*>(c->Clone());
298 fShapeCorr->SetDirectory(0);
301 //________________________________________________________________________
303 AliBasedNdetaTask::InitializeCentBins()
305 if (fListOfCentralities->GetEntries() > 0) return;
307 // Automatically add 'all' centrality bin if nothing has been defined.
308 AddCentralityBin(0, 0, 0);
309 if (fCentAxis && fCentAxis->GetNbins() > 0 && fCentAxis->GetXbins()) {
310 const TArrayD* bins = fCentAxis->GetXbins();
311 Int_t nbin = fCentAxis->GetNbins();
312 for (Int_t i = 0; i < nbin; i++)
313 AddCentralityBin(i+1, Short_t((*bins)[i]), Short_t((*bins)[i+1]));
317 //________________________________________________________________________
319 AliBasedNdetaTask::UserCreateOutputObjects()
322 // Create output objects.
324 // This is called once per slave process
326 DGUARD(fDebug,1,"Create user ouput object");
328 fSums->SetName(Form("%s_sums", GetName()));
331 InitializeCentBins();
332 if (fCentAxis) fSums->Add(fCentAxis);
335 // Centrality histogram
336 fCent = new TH1D("cent", "Centrality", 100, 0, 100);
337 fCent->SetDirectory(0);
341 // Loop over centrality bins
342 TIter next(fListOfCentralities);
343 CentralityBin* bin = 0;
344 while ((bin = static_cast<CentralityBin*>(next())))
345 bin->CreateOutputObjects(fSums);
348 // Post data for ALL output slots >0 here, to get at least an empty
353 //____________________________________________________________________
355 AliBasedNdetaTask::UserExec(Option_t *)
358 // Process a single event
364 DGUARD(fDebug,1,"Analyse the AOD event");
365 AliAODEvent* aod = AliForwardUtil::GetAODEvent(this);
369 TObject* obj = aod->FindListObject("Forward");
371 AliWarning("No forward object found");
374 AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
376 // Fill centrality histogram
377 Float_t cent = forward->GetCentrality();
380 // Get the histogram(s)
381 TH2D* data = GetHistogram(aod, false);
382 TH2D* dataMC = GetHistogram(aod, true);
384 Bool_t isZero = ((fNormalizationScheme & kZeroBin) &&
385 !forward->IsTriggerBits(AliAODForwardMult::kNClusterGt0));
388 // Loop over centrality bins
389 CentralityBin* allBin =
390 static_cast<CentralityBin*>(fListOfCentralities->At(0));
391 allBin->ProcessEvent(forward, fTriggerMask, isZero,
392 fVtxMin, fVtxMax, data, dataMC);
394 // Find this centrality bin
395 if (fCentAxis && fCentAxis->GetNbins() > 0) {
396 Int_t icent = fCentAxis->FindBin(cent);
397 CentralityBin* thisBin = 0;
398 if (icent >= 1 && icent <= fCentAxis->GetNbins())
399 thisBin = static_cast<CentralityBin*>(fListOfCentralities->At(icent));
401 thisBin->ProcessEvent(forward, fTriggerMask, isZero, fVtxMin, fVtxMax,
405 // Here, we get the update
407 UShort_t sNN = forward->GetSNN();
408 fSNNString = new TNamed("sNN", "");
409 fSNNString->SetTitle(AliForwardUtil::CenterOfMassEnergyString(sNN));
410 fSNNString->SetUniqueID(sNN);
411 fSums->Add(fSNNString);
413 UShort_t sys = forward->GetSystem();
414 fSysString = new TNamed("sys", "");
415 fSysString->SetTitle(AliForwardUtil::CollisionSystemString(sys));
416 fSysString->SetUniqueID(sys);
417 fSums->Add(fSysString);
419 fSums->Add(fSchemeString);
420 fSums->Add(fTriggerString);
428 //________________________________________________________________________
430 AliBasedNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker,
431 const char* title, const char* ytitle)
434 // Set histogram graphical options, etc.
437 // h Histogram to modify
438 // colour Marker color
439 // marker Marker style
440 // title Title of histogram
441 // ytitle Title on y-axis.
444 h->SetMarkerColor(colour);
445 h->SetMarkerStyle(marker);
446 h->SetMarkerSize(marker == 29 || marker == 30 ? 1.2 : 1);
448 h->SetYTitle(ytitle);
449 h->GetXaxis()->SetTitleFont(132);
450 h->GetXaxis()->SetLabelFont(132);
451 h->GetXaxis()->SetNdivisions(10);
452 h->GetYaxis()->SetTitleFont(132);
453 h->GetYaxis()->SetLabelFont(132);
454 h->GetYaxis()->SetNdivisions(10);
455 h->GetYaxis()->SetDecimals();
459 //________________________________________________________________________
461 AliBasedNdetaTask::ScaleToCoverage(TH2D* copy, const TH1D* norm)
463 // Normalize to the acceptance -
464 // dndeta->Divide(accNorm);
465 for (Int_t i = 1; i <= copy->GetNbinsX(); i++) {
466 Double_t a = norm->GetBinContent(i);
467 for (Int_t j = 1; j <= copy->GetNbinsY(); j++) {
469 copy->SetBinContent(i,j,0);
470 copy->SetBinError(i,j,0);
473 Double_t c = copy->GetBinContent(i, j);
474 Double_t e = copy->GetBinError(i, j);
475 copy->SetBinContent(i, j, c / a);
476 copy->SetBinError(i, j, e / a);
481 //________________________________________________________________________
483 AliBasedNdetaTask::ProjectX(const TH2D* h,
492 // Project onto the X axis
497 // firstbin First bin to use
498 // lastbin Last bin to use
499 // error Whether to calculate errors
502 // Newly created histogram or null
506 return h->ProjectionX(name, firstbin, lastbin, (error ? "e" : ""));
508 TAxis* xaxis = h->GetXaxis();
509 TAxis* yaxis = h->GetYaxis();
510 TH1D* ret = new TH1D(name, h->GetTitle(), xaxis->GetNbins(),
511 xaxis->GetXmin(), xaxis->GetXmax());
512 static_cast<const TAttLine*>(h)->Copy(*ret);
513 static_cast<const TAttFill*>(h)->Copy(*ret);
514 static_cast<const TAttMarker*>(h)->Copy(*ret);
515 ret->GetXaxis()->ImportAttributes(xaxis);
517 Int_t first = firstbin;
518 Int_t last = lastbin;
519 if (first < 0) first = 0;
520 else if (first >= yaxis->GetNbins()+1) first = yaxis->GetNbins();
521 if (last < 0) last = yaxis->GetNbins();
522 else if (last > yaxis->GetNbins()+1) last = yaxis->GetNbins();
523 if (last-first < 0) {
524 AliWarningGeneral("AliBasedNdetaTask",
525 Form("Nothing to project [%d,%d]", first, last));
531 // AliInfo(Form("Projecting bins [%d,%d] of %s", first, last, h->GetName()));
532 Int_t ybins = (last-first+1);
533 for (Int_t xbin = 0; xbin <= xaxis->GetNbins()+1; xbin++) {
534 Double_t content = 0;
539 for (Int_t ybin = first; ybin <= last; ybin++) {
540 Double_t c1 = h->GetCellContent(xbin, ybin);
541 Double_t e1 = h->GetCellError(xbin, ybin);
544 if (c1 < 1e-12) continue;
554 if(content > 0 && nbins > 0) {
555 Double_t factor = (corr ? Double_t(ybins) / nbins : 1);
557 AliWarningGeneral(ret->GetName(),
558 Form("factor @ %d is %d/%d -> %f",
559 xbin, ybins, nbins, factor));
562 // Calculate weighted average
563 ret->SetBinContent(xbin, content * factor);
564 ret->SetBinError(xbin, factor * TMath::Sqrt(error2));
567 ret->SetBinContent(xbin, factor * content);
574 //________________________________________________________________________
576 AliBasedNdetaTask::Terminate(Option_t *)
579 // Called at end of event processing..
581 // This is called once in the master
586 // Draw result to screen, or perform fitting, normalizations Called
587 // once at the end of the query
588 DGUARD(fDebug,1,"Process final merged results");
590 fSums = dynamic_cast<TList*> (GetOutputData(1));
592 AliError("Could not retrieve TList fSums");
597 fOutput->SetName(Form("%s_result", GetName()));
600 fSNNString = static_cast<TNamed*>(fSums->FindObject("sNN"));
601 fSysString = static_cast<TNamed*>(fSums->FindObject("sys"));
602 fCentAxis = static_cast<TAxis*>(fSums->FindObject("centAxis"));
603 fSchemeString = static_cast<TNamed*>(fSums->FindObject("scheme"));
604 fTriggerString = static_cast<TNamed*>(fSums->FindObject("trigger"));
606 if(fSysString && fSNNString &&
607 fSysString->GetUniqueID() == AliForwardUtil::kPP)
608 LoadNormalizationData(fSysString->GetUniqueID(),
609 fSNNString->GetUniqueID());
611 InitializeCentBins();
613 // Print before we loop
616 // Loop over centrality bins
617 TIter next(fListOfCentralities);
618 CentralityBin* bin = 0;
619 gStyle->SetPalette(1);
620 THStack* dndetaStack = new THStack("dndeta", "dN/d#eta");
621 THStack* dndetaStackRebin = new THStack(Form("dndeta_rebin%02d", fRebin),
623 THStack* dndetaMCStack = new THStack("dndetaMC", "dN_{ch}/d#eta");
624 THStack* dndetaMCStackRebin = new THStack(Form("dndetaMC_rebin%02d", fRebin),
628 TList* truthlist = 0;
630 if (fFinalMCCorrFile.Contains(".root")) {
631 TFile* ftest = TFile::Open(fFinalMCCorrFile.Data());
633 mclist = dynamic_cast<TList*>(ftest->Get(Form("%sResults",GetName())));
634 truthlist = dynamic_cast<TList*>(ftest->Get("MCTruthResults"));
637 AliWarning("MC analysis file invalid - no final MC correction possible");
639 Int_t style = GetMarker();
640 Int_t color = GetColor();
642 AliInfo(Form("Marker style=%d, color=%d", style, color));
643 while ((bin = static_cast<CentralityBin*>(next()))) {
645 bin->End(fSums, fOutput, fNormalizationScheme, fShapeCorr, fTriggerEff,
646 fSymmetrice, fRebin, fUseROOTProj, fCorrEmpty, fCutEdges,
647 fTriggerMask, style, color, mclist, truthlist);
648 if (fCentAxis && bin->IsAllBin()) continue;
649 TH1* dndeta = bin->GetResult(0, false, "");
650 TH1* dndetaSym = bin->GetResult(0, true, "");
651 TH1* dndetaMC = bin->GetResult(0, false, "MC");
652 TH1* dndetaMCSym = bin->GetResult(0, true, "MC");
653 if (dndeta) dndetaStack->Add(dndeta);
654 if (dndetaSym) dndetaStack->Add(dndetaSym);
655 if (dndetaMC) dndetaMCStack->Add(dndetaMC);
656 if (dndetaMCSym) dndetaMCStack->Add(dndetaMCSym);
658 dndeta = bin->GetResult(fRebin, false, "");
659 dndetaSym = bin->GetResult(fRebin, true, "");
660 dndetaMC = bin->GetResult(fRebin, false, "MC");
661 dndetaMCSym = bin->GetResult(fRebin, true, "MC");
662 if (dndeta) dndetaStackRebin->Add(dndeta);
663 if (dndetaSym) dndetaStackRebin->Add(dndetaSym);
664 if (dndetaMC) dndetaMCStackRebin->Add(dndetaMC);
665 if (dndetaMCSym) dndetaMCStackRebin->Add(dndetaMCSym);
669 fOutput->Add(dndetaStack);
671 // If available output rebinned stack
672 if (!dndetaStackRebin->GetHists() ||
673 dndetaStackRebin->GetHists()->GetEntries() <= 0) {
674 AliWarning("No rebinned histograms found");
675 delete dndetaStackRebin;
676 dndetaStackRebin = 0;
678 if (dndetaStackRebin) fOutput->Add(dndetaStackRebin);
680 // If available, output track-ref stack
681 if (!dndetaMCStack->GetHists() ||
682 dndetaMCStack->GetHists()->GetEntries() <= 0) {
683 AliWarning("No MC histograms found");
684 delete dndetaMCStack;
687 if (dndetaMCStack) fOutput->Add(dndetaMCStack);
689 // If available, output rebinned track-ref stack
690 if (!dndetaMCStackRebin->GetHists() ||
691 dndetaMCStackRebin->GetHists()->GetEntries() <= 0) {
692 AliWarning("No rebinned MC histograms found");
693 delete dndetaMCStackRebin;
694 dndetaMCStackRebin = 0;
696 if (dndetaMCStackRebin) fOutput->Add(dndetaMCStackRebin);
698 // Output collision energy string
699 if (fSNNString) fOutput->Add(fSNNString->Clone());
701 // Output collision system string
702 if (fSysString) fOutput->Add(fSysString->Clone());
704 // Output centrality axis
705 if (fCentAxis) fOutput->Add(fCentAxis);
707 // Output trigger string
708 if (fTriggerString) fOutput->Add(fTriggerString->Clone());
710 // Normalization string
711 if (fSchemeString) fOutput->Add(fSchemeString->Clone());
713 // Output vertex axis
714 TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax);
715 vtxAxis->SetName("vtxAxis");
716 vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
717 fOutput->Add(vtxAxis);
719 // Output trigger efficiency and shape correction
720 fOutput->Add(new TParameter<Double_t>("triggerEff", fTriggerEff));
721 if (fShapeCorr) fOutput->Add(fShapeCorr);
723 TNamed* options = new TNamed("options","");
725 str.Append(Form("Edges %scut, ", fCutEdges ? "" : "not "));
726 str.Append(Form("Empty bins %scorrected for, ", fCorrEmpty ? "" : "not "));
727 str.Append(Form("TH2::ProjectionX %sused", fUseROOTProj ? "" : "not "));
728 options->SetTitle(str);
729 fOutput->Add(options);
731 PostData(2, fOutput);
733 //________________________________________________________________________
735 AliBasedNdetaTask::LoadNormalizationData(UShort_t sys, UShort_t energy)
737 // Load the normalisation data for dN/deta for pp INEL, INEL>0, and NSD
738 DGUARD(fDebug,1,"Load the normalization data for sys=%d, energy=%d",
742 if(energy == 7000) snn.Form("7000");
743 if(energy == 2750) snn.Form("2750");
745 if(fShapeCorr && (fTriggerEff != 1)) {
746 AliInfo("Objects already set for normalization - no action taken");
750 TFile* fin = TFile::Open(Form("$ALICE_ROOT/PWGLF/FORWARD/corrections/"
751 "Normalization/normalizationHists_%s_%s.root",
752 type.Data(),snn.Data()));
754 AliWarning(Form("no file for normalization of %d/%d", sys, energy));
759 if ((fNormalizationScheme & kShape) && !fShapeCorr) {
760 TString trigName("All");
761 if (fTriggerMask == AliAODForwardMult::kInel ||
762 fTriggerMask == AliAODForwardMult::kNClusterGt0)
764 else if (fTriggerMask == AliAODForwardMult::kNSD)
766 else if (fTriggerMask == AliAODForwardMult::kInelGt0)
767 trigName = "InelGt0";
769 AliWarning(Form("Normalization for trigger %s not known, using all",
770 AliAODForwardMult::GetTriggerString(fTriggerMask)));
773 TString shapeCorName(Form("h%sNormalization", trigName.Data()));
774 TH2F* shapeCor = dynamic_cast<TH2F*>(fin->Get(shapeCorName));
775 if (shapeCor) SetShapeCorrection(shapeCor);
777 AliWarning(Form("No shape correction found for %s", trigName.Data()));
781 // Trigger efficiency
782 TString effName(Form("%sTriggerEff",
783 fTriggerMask == AliAODForwardMult::kInel ? "inel" :
784 fTriggerMask == AliAODForwardMult::kNSD ? "nsd" :
785 fTriggerMask == AliAODForwardMult::kInelGt0 ?
787 TParameter<float>* eff = 0;
788 if (fNormalizationScheme & kTriggerEfficiency)
789 eff = static_cast<TParameter<float>*>(fin->Get(effName));
790 Double_t trigEff = eff ? eff->GetVal() : 1;
791 if (fTriggerEff != 1) SetTriggerEff(trigEff);
792 if (fTriggerEff < 0) fTriggerEff = 1;
795 // Rescale the shape correction by the trigger efficiency
797 AliWarning(Form("Rescaling shape correction by trigger efficiency: "
798 "1/E_X=1/%f", trigEff));
799 fShapeCorr->Scale(1. / trigEff);
803 if (fShapeCorr && fTriggerEff) AliInfo("Loaded objects for normalization.");
807 //________________________________________________________________________
809 AliBasedNdetaTask::Print(Option_t*) const
814 std::cout << this->ClassName() << ": " << this->GetName() << "\n"
816 << " Trigger: " << (fTriggerString ?
817 fTriggerString->GetTitle() :
819 << " Vertex range: [" << fVtxMin << ":" << fVtxMax << "]\n"
820 << " Rebin factor: " << fRebin << "\n"
821 << " Cut edges: " << fCutEdges << "\n"
822 << " Symmertrice: " << fSymmetrice << "\n"
823 << " Use TH2::ProjectionX: " << fUseROOTProj << "\n"
824 << " Correct for empty: " << fCorrEmpty << "\n"
825 << " Normalization scheme: " << (fSchemeString ?
826 fSchemeString->GetTitle() :
828 << " Trigger efficiency: " << fTriggerEff << "\n"
829 << " Shape correction: " << (fShapeCorr ?
830 fShapeCorr->GetName() :
832 << " sqrt(s_NN): " << (fSNNString ?
833 fSNNString->GetTitle() :
835 << " Collision system: " << (fSysString ?
836 fSysString->GetTitle() :
838 << " Centrality bins: " << (fCentAxis ? "" : "none");
840 Int_t nBins = fCentAxis->GetNbins();
841 const Double_t* bins = fCentAxis->GetXbins()->GetArray();
842 for (Int_t i = 0; i <= nBins; i++)
843 std::cout << (i==0 ? "" : "-") << bins[i];
845 std::cout << std::noboolalpha << std::endl;
849 //________________________________________________________________________
851 AliBasedNdetaTask::Rebin(const TH1D* h, Int_t rebin, Bool_t cutEdges)
854 // Make a copy of the input histogram and rebin that histogram
857 // h Histogram to rebin
860 // New (rebinned) histogram
862 if (rebin <= 1) return 0;
864 Int_t nBins = h->GetNbinsX();
865 if(nBins % rebin != 0) {
866 AliWarningGeneral("AliBasedNdetaTask",
867 Form("Rebin factor %d is not a devisor of current number "
868 "of bins %d in the histogram %s",
869 rebin, nBins, h->GetName()));
874 TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d",
875 h->GetName(), rebin)));
877 tmp->SetDirectory(0);
879 // The new number of bins
880 Int_t nBinsNew = nBins / rebin;
881 for(Int_t i = 1;i<= nBinsNew; i++) {
882 Double_t content = 0;
886 for(Int_t j = 1; j<=rebin;j++) {
887 Int_t bin = (i-1)*rebin + j;
888 Double_t c = h->GetBinContent(bin);
889 if (c <= 0) continue;
892 if (h->GetBinContent(bin+1)<=0 ||
893 h->GetBinContent(bin-1)<=0) {
895 AliWarningGeneral("AliBasedNdetaTask",
896 Form("removing bin %d=%f of %s (%d=%f,%d=%f)",
897 bin, c, h->GetName(),
898 bin+1, h->GetBinContent(bin+1),
899 bin-1, h->GetBinContent(bin-1)));
904 Double_t e = h->GetBinError(bin);
905 Double_t w = 1 / (e*e); // 1/c/c
912 if(content > 0 && nbins > 0) {
913 tmp->SetBinContent(i, wsum / sumw);
914 tmp->SetBinError(i,1./TMath::Sqrt(sumw));
921 //__________________________________________________________________
923 AliBasedNdetaTask::Symmetrice(const TH1* h)
926 // Make an extension of @a h to make it symmetric about 0
929 // h Histogram to symmertrice
932 // Symmetric extension of @a h
934 Int_t nBins = h->GetNbinsX();
935 TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
936 s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
938 s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
939 s->SetMarkerColor(h->GetMarkerColor());
940 s->SetMarkerSize(h->GetMarkerSize());
941 s->SetMarkerStyle(FlipHollowStyle(h->GetMarkerStyle()));
942 s->SetFillColor(h->GetFillColor());
943 s->SetFillStyle(h->GetFillStyle());
946 // Find the first and last bin with data
947 Int_t first = nBins+1;
949 for (Int_t i = 1; i <= nBins; i++) {
950 if (h->GetBinContent(i) <= 0) continue;
951 first = TMath::Min(first, i);
952 last = TMath::Max(last, i);
955 Double_t xfirst = h->GetBinCenter(first-1);
956 Int_t f1 = h->GetXaxis()->FindBin(-xfirst);
957 Int_t l2 = s->GetXaxis()->FindBin(xfirst);
958 for (Int_t i = f1, j=l2; i <= last; i++,j--) {
959 s->SetBinContent(j, h->GetBinContent(i));
960 s->SetBinError(j, h->GetBinError(i));
962 // Fill in overlap bin
963 s->SetBinContent(l2+1, h->GetBinContent(first));
964 s->SetBinError(l2+1, h->GetBinError(first));
968 //__________________________________________________________________
970 AliBasedNdetaTask::GetMarkerStyle(UShort_t bits)
972 Int_t base = bits & (0xFE);
973 Bool_t hollow = bits & kHollow;
975 case kCircle: return (hollow ? 24 : 20);
976 case kSquare: return (hollow ? 25 : 21);
977 case kUpTriangle: return (hollow ? 26 : 22);
978 case kDownTriangle: return (hollow ? 32 : 23);
979 case kDiamond: return (hollow ? 27 : 33);
980 case kCross: return (hollow ? 28 : 34);
981 case kStar: return (hollow ? 30 : 29);
985 //__________________________________________________________________
987 AliBasedNdetaTask::GetMarkerBits(Int_t style)
991 case 24: case 25: case 26: case 27: case 28: case 30: case 32:
992 bits |= kHollow; break;
995 case 20: case 24: bits |= kCircle; break;
996 case 21: case 25: bits |= kSquare; break;
997 case 22: case 26: bits |= kUpTriangle; break;
998 case 23: case 32: bits |= kDownTriangle; break;
999 case 27: case 33: bits |= kDiamond; break;
1000 case 28: case 34: bits |= kCross; break;
1001 case 29: case 30: bits |= kStar; break;
1005 //__________________________________________________________________
1007 AliBasedNdetaTask::FlipHollowStyle(Int_t style)
1009 UShort_t bits = GetMarkerBits(style);
1010 Int_t ret = GetMarkerStyle(bits ^ kHollow);
1014 //====================================================================
1016 AliBasedNdetaTask::Sum::Init(TList* list, const TH2D* data, Int_t col)
1018 DGUARD(fDebug,1,"Initializing sums with %s", data->GetName());
1019 TString n(GetHistName(0));
1020 TString n0(GetHistName(1));
1021 const char* postfix = GetTitle();
1023 fSum = static_cast<TH2D*>(data->Clone(n));
1024 if (postfix) fSum->SetTitle(Form("%s (%s)", data->GetTitle(), postfix));
1025 fSum->SetDirectory(0);
1026 fSum->SetMarkerColor(col);
1027 fSum->SetMarkerStyle(GetMarkerStyle(kCircle|kSolid));
1031 fSum0 = static_cast<TH2D*>(data->Clone(n0));
1033 fSum0->SetTitle(Form("%s 0-bin (%s)", data->GetTitle(), postfix));
1035 fSum0->SetTitle(Form("%s 0-bin", data->GetTitle()));
1036 fSum0->SetDirectory(0);
1037 fSum0->SetMarkerColor(col);
1038 fSum0->SetMarkerStyle(GetMarkerStyle(kCross|kHollow));
1042 fEvents = new TH1I(GetHistName(2), "Event types", 2, -.5, 1.5);
1043 fEvents->SetDirectory(0);
1044 fEvents->GetXaxis()->SetBinLabel(1, "Non-zero");
1045 fEvents->GetXaxis()->SetBinLabel(2, "Zero");
1049 //____________________________________________________________________
1051 AliBasedNdetaTask::Sum::GetHistName(const char* name, Int_t what, const char* post)
1054 if (what == 1) n.Append("0");
1055 else if (what == 2) n.Append("Events");
1056 if (post && post[0] != '\0') n.Append(post);
1060 //____________________________________________________________________
1062 AliBasedNdetaTask::Sum::GetHistName(Int_t what) const
1064 return GetHistName(GetName(), what, GetTitle());
1067 //____________________________________________________________________
1069 AliBasedNdetaTask::Sum::Add(const TH2D* data, Bool_t isZero)
1071 DGUARD(fDebug,2,"Adding %s to sums", data->GetName());
1072 if (isZero) fSum0->Add(data);
1073 else fSum->Add(data);
1074 fEvents->Fill(isZero ? 1 : 0);
1077 //____________________________________________________________________
1079 AliBasedNdetaTask::Sum::CalcSum(TList* output,
1085 Bool_t corrEmpty) const
1087 DGUARD(fDebug,2,"Calculating final summed histogram %s", fSum->GetName());
1089 TH2D* ret = static_cast<TH2D*>(fSum->Clone(fSum->GetName()));
1090 ret->SetDirectory(0);
1092 Int_t n = Int_t(fEvents->GetBinContent(1));
1093 Int_t n0 = Int_t(fEvents->GetBinContent(2));
1095 AliInfoF("Adding histograms %s and %s with weights %f and %f resp.",
1096 fSum0->GetName(), fSum->GetName(), 1./epsilon, 1./epsilon0);
1097 DMSG(fDebug,2, "Adding histograms %s and %s with weights %f and %f resp.",
1098 fSum0->GetName(), fSum->GetName(), 1./epsilon, 1./epsilon0);
1099 // Generate merged histogram
1100 ret->Add(fSum0, fSum, 1. / epsilon0, 1. / epsilon);
1101 ntotal = n / epsilon + n0 / epsilon0;
1103 TList* out = new TList;
1105 const char* postfix = GetTitle();
1106 if (!postfix) postfix = "";
1107 out->SetName(Form("partial%s", postfix));
1110 // Now make copies, normalize them, and store in output list
1111 TH2D* sumCopy = static_cast<TH2D*>(fSum->Clone("sum"));
1112 TH2D* sum0Copy = static_cast<TH2D*>(fSum0->Clone("sum0"));
1113 TH2D* retCopy = static_cast<TH2D*>(ret->Clone("sumAll"));
1114 sumCopy->SetMarkerStyle(FlipHollowStyle(marker));
1115 sumCopy->SetDirectory(0);
1116 sum0Copy->SetMarkerStyle(GetMarkerStyle(GetMarkerBits(marker)+4));
1117 sum0Copy->SetDirectory(0);
1118 retCopy->SetMarkerStyle(marker);
1119 retCopy->SetDirectory(0);
1121 TH1D* norm = ProjectX(fSum, "norm", 0, 0, rootProj, corrEmpty, false);
1122 TH1D* norm0 = ProjectX(fSum0, "norm0", 0, 0, rootProj, corrEmpty, false);
1123 TH1D* normAll = ProjectX(ret, "normAll", 0, 0, rootProj, corrEmpty, false);
1124 norm->SetDirectory(0);
1125 norm0->SetDirectory(0);
1126 normAll->SetDirectory(0);
1128 ScaleToCoverage(sumCopy, norm);
1129 ScaleToCoverage(sum0Copy, norm0);
1130 ScaleToCoverage(retCopy, normAll);
1132 Int_t nY = fSum->GetNbinsY();
1133 TH1D* sumCopyPx = ProjectX(sumCopy, "average", 1, nY,rootProj,corrEmpty);
1134 TH1D* sum0CopyPx = ProjectX(sum0Copy, "average0", 1, nY,rootProj,corrEmpty);
1135 TH1D* retCopyPx = ProjectX(retCopy, "averageAll", 1, nY,rootProj,corrEmpty);
1136 sumCopyPx->SetDirectory(0);
1137 sum0CopyPx->SetDirectory(0);
1138 retCopyPx->SetDirectory(0);
1140 // Scale our 1D histograms
1141 sumCopyPx->Scale(1., "width");
1142 sum0CopyPx->Scale(1., "width");
1143 retCopyPx->Scale(1., "width");
1145 AliInfo(Form("Maximum %f,%f changed to %f", sumCopyPx->GetMaximum(),
1146 sum0CopyPx->GetMaximum(), retCopyPx->GetMaximum()));
1148 // Scale the normalization - they should be 1 at the maximum
1149 norm->Scale(n > 0 ? 1. / n : 1);
1150 norm0->Scale(n0 > 0 ? 1. / n0 : 1);
1151 normAll->Scale(ntotal > 0 ? 1. / ntotal : 1);
1156 out->Add(sumCopyPx);
1157 out->Add(sum0CopyPx);
1158 out->Add(retCopyPx);
1163 AliInfoF("Returning (1/%f * %s + 1/%f * %s), "
1164 "1/%f * %d + 1/%f * %d = %d",
1165 epsilon0, fSum0->GetName(), epsilon, fSum->GetName(),
1166 epsilon0, n0, epsilon, n, int(ntotal));
1168 for (Int_t i = 1; i <= ret->GetNbinsX(); i++) {
1169 Double_t nc = sum->GetBinContent(i, 0);
1170 Double_t nc0 = sum0->GetBinContent(i, 0);
1171 ret->SetBinContent(i, 0, nc + nc0); // Just count events
1178 //====================================================================
1179 AliBasedNdetaTask::CentralityBin::CentralityBin()
1188 fDoFinalMCCorrection(false),
1194 DGUARD(0,0,"Centrality bin default construction");
1196 //____________________________________________________________________
1197 AliBasedNdetaTask::CentralityBin::CentralityBin(const char* name,
1198 Short_t low, Short_t high)
1207 fDoFinalMCCorrection(false),
1214 // name Name used for histograms (e.g., Forward)
1215 // low Lower centrality cut in percent
1216 // high Upper centrality cut in percent
1218 DGUARD(0,0,"Named Centrality bin construction: %s [%d,%d]", name);
1219 if (low <= 0 && high <= 0) {
1222 SetTitle("All centralities");
1227 SetTitle(Form("Centrality bin from %3d%% to %3d%%", low, high));
1230 //____________________________________________________________________
1231 AliBasedNdetaTask::CentralityBin::CentralityBin(const CentralityBin& o)
1237 fTriggers(o.fTriggers),
1240 fDoFinalMCCorrection(o.fDoFinalMCCorrection),
1247 // other Object to copy from
1249 DGUARD(0,0,"Copy Centrality bin construction");
1251 //____________________________________________________________________
1252 AliBasedNdetaTask::CentralityBin::~CentralityBin()
1257 DGUARD(fDebug,3,"Centrality bin desctruction");
1258 if (fSums) fSums->Delete();
1259 if (fOutput) fOutput->Delete();
1262 //____________________________________________________________________
1263 AliBasedNdetaTask::CentralityBin&
1264 AliBasedNdetaTask::CentralityBin::operator=(const CentralityBin& o)
1267 // Assignment operator
1270 // other Object to assign from
1273 // Reference to this
1275 DGUARD(fDebug,3,"Centrality bin assignment");
1276 if (&o == this) return *this;
1277 SetName(o.GetName());
1278 SetTitle(o.GetTitle());
1280 fOutput = o.fOutput;
1283 fTriggers = o.fTriggers;
1286 fDoFinalMCCorrection = o.fDoFinalMCCorrection;
1290 //____________________________________________________________________
1292 AliBasedNdetaTask::CentralityBin::GetColor(Int_t fallback) const
1294 if (IsAllBin()) return fallback;
1295 Float_t fc = (fLow+double(fHigh-fLow)/2) / 100;
1296 Int_t nCol = gStyle->GetNumberOfColors();
1297 Int_t icol = TMath::Min(nCol-1,int(fc * nCol + .5));
1298 Int_t col = gStyle->GetColorPalette(icol);
1301 //____________________________________________________________________
1303 AliBasedNdetaTask::CentralityBin::GetListName() const
1306 // Get the list name
1311 if (IsAllBin()) return "all";
1312 return Form("cent%03d_%03d", fLow, fHigh);
1314 //____________________________________________________________________
1316 AliBasedNdetaTask::CentralityBin::CreateOutputObjects(TList* dir)
1319 // Create output objects
1324 DGUARD(fDebug,1,"Create centrality bin output objects");
1326 fSums->SetName(GetListName());
1330 fTriggers = AliAODForwardMult::MakeTriggerHistogram("triggers");
1331 fTriggers->SetDirectory(0);
1332 fSums->Add(fTriggers);
1334 //____________________________________________________________________
1336 AliBasedNdetaTask::CentralityBin::SetDebugLevel(Int_t lvl)
1339 if (fSum) fSum->fDebug = lvl;
1340 if (fSumMC) fSumMC->fDebug = lvl;
1343 //____________________________________________________________________
1345 AliBasedNdetaTask::CentralityBin::ReadSum(TList* list, bool mc)
1347 const char* post = (mc ? "MC" : "");
1348 TString sn = Sum::GetHistName(GetName(),0,post);
1349 TString sn0 = Sum::GetHistName(GetName(),1,post);
1350 TString ev = Sum::GetHistName(GetName(),2,post);
1351 TH2D* sum = static_cast<TH2D*>(list->FindObject(sn));
1352 TH2D* sum0 = static_cast<TH2D*>(list->FindObject(sn0));
1353 TH1I* events = static_cast<TH1I*>(list->FindObject(ev));
1354 if (!sum || !sum0 || !events) {
1355 AliWarningF("Failed to find one or more histograms: "
1356 "%s (%p) %s (%p) %s (%p)",
1362 Sum* ret = new Sum(GetName(), post);
1365 ret->fEvents = events;
1366 ret->fDebug = fDebug;
1367 if (mc) fSumMC = ret;
1373 //____________________________________________________________________
1375 AliBasedNdetaTask::CentralityBin::CreateSums(const TH2D* data, const TH2D* mc)
1378 // Create sum histogram
1381 // data Data histogram to clone
1382 // mc (optional) MC histogram to clone
1384 DGUARD(fDebug,1,"Create centrality bin sums from %s",
1385 data ? data->GetName() : "(null)");
1387 fSum = new Sum(GetName(),"");
1388 fSum->Init(fSums, data, GetColor());
1389 fSum->fDebug = fDebug;
1392 // If no MC data is given, then do not create MC sum histogram
1395 fSumMC = new Sum(GetName(), "MC");
1396 fSumMC->Init(fSums, mc, GetColor());
1397 fSumMC->fDebug = fDebug;
1400 //____________________________________________________________________
1402 AliBasedNdetaTask::CentralityBin::CheckEvent(const AliAODForwardMult* forward,
1404 Double_t vzMin, Double_t vzMax)
1407 // Check the trigger, vertex, and centrality
1413 // true if the event is to be used
1415 if (!forward) return false;
1417 DGUARD(fDebug,2,"Check the event");
1418 // We do not check for centrality here - it's already done
1419 return forward->CheckEvent(triggerMask, vzMin, vzMax, 0, 0, fTriggers);
1423 //____________________________________________________________________
1425 AliBasedNdetaTask::CentralityBin::ProcessEvent(const AliAODForwardMult* forward,
1426 Int_t triggerMask, Bool_t isZero,
1427 Double_t vzMin, Double_t vzMax,
1428 const TH2D* data, const TH2D* mc)
1434 // forward Forward data (for trigger, vertex, & centrality)
1435 // triggerMask Trigger mask
1436 // vzMin Minimum IP z coordinate
1437 // vzMax Maximum IP z coordinate
1438 // data Data histogram
1441 DGUARD(fDebug,1,"Process one event for %s a given centrality bin",
1442 data ? data->GetName() : "(null)");
1443 if (!CheckEvent(forward, triggerMask, vzMin, vzMax)) return;
1445 if (!fSum) CreateSums(data, mc);
1447 fSum->Add(data, isZero);
1448 if (mc) fSumMC->Add(mc, isZero);
1451 //________________________________________________________________________
1453 AliBasedNdetaTask::CentralityBin::Normalization(const TH1I& t,
1456 Double_t& ntotal) const
1459 // Calculate normalization
1462 // t Trigger histogram
1463 // scheme Normaliztion scheme
1465 // ntotal On return, contains the number of events.
1467 DGUARD(fDebug,1,"Normalize centrality bin with %s", t.GetName());
1468 Double_t nAll = t.GetBinContent(AliAODForwardMult::kBinAll);
1469 Double_t nB = t.GetBinContent(AliAODForwardMult::kBinB);
1470 Double_t nA = t.GetBinContent(AliAODForwardMult::kBinA);
1471 Double_t nC = t.GetBinContent(AliAODForwardMult::kBinC);
1472 Double_t nE = t.GetBinContent(AliAODForwardMult::kBinE);
1473 Double_t nOffline = t.GetBinContent(AliAODForwardMult::kBinOffline);
1474 Double_t nTriggered = t.GetBinContent(AliAODForwardMult::kWithTrigger);
1475 Double_t nWithVertex = t.GetBinContent(AliAODForwardMult::kWithVertex);
1476 Double_t nAccepted = ntotal; // t.GetBinContent(AliAODForwardMult::kAccepted);
1479 if (nTriggered <= 0.1) {
1480 AliError("Number of triggered events <= 0");
1483 if (nWithVertex <= 0.1) {
1484 AliError("Number of events with vertex <= 0");
1488 Double_t vtxEff = nWithVertex / nTriggered;
1489 Double_t scaler = 1;
1490 Double_t beta = nA + nC - 2*nE;
1492 if (scheme & kEventLevel && !(scheme & kZeroBin)) {
1493 ntotal = nAccepted / vtxEff;
1495 AliInfo(Form("Calculating event normalisation as\n"
1496 " N = N_A * N_T / N_V = %d * %d / %d = %f (%f)",
1497 Int_t(nAccepted), Int_t(nTriggered), Int_t(nWithVertex),
1500 if (scheme & kBackground) {
1502 // s = --------- = ------------- = ------------
1503 // 1 - beta 1 - beta E_V 1 - beta N_V
1504 // --- ---- -------- ---- ---
1505 // E_V N_V N_V N_V N_T
1513 ntotal -= nAccepted * beta / nWithVertex;
1514 // This one is direct and correct.
1515 // scaler = 1. / (1. / vtxEff - beta / nWithVertex);
1516 // A simpler expresion
1517 scaler /= (1 - beta / nTriggered); // 0.831631 -> 0.780689
1518 AliInfo(Form("Calculating event normalisation as\n"
1519 " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n"
1520 " N = N - N_A * beta / N_V = %f - %d * %d / %d = %f (%f)",
1521 Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta),
1522 nAccepted / vtxEff, Int_t(nAccepted), Int_t(beta),
1523 Int_t(nWithVertex), ntotal, scaler));
1526 if (scheme & kZeroBin) {
1529 // N = N_A + 1/E_X * N_A / N_V (N_T - N_V - beta)
1530 // = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V))
1532 // s = N_A/N = 1 / (1 + 1/E_X (N_T/N_V - 1 - beta / N_V))
1533 // = N_V / (N_V + 1/E_X (N_T - N_V - beta))
1535 if (!(scheme & kBackground)) beta = 0;
1536 ntotal = nAccepted * (1 + 1/trigEff * (nTriggered / nWithVertex - 1
1537 - beta / nWithVertex));
1538 scaler = nWithVertex / (nWithVertex +
1539 1/trigEff * (nTriggered-nWithVertex-beta));
1540 AliInfo(Form("Calculating event normalisation as\n"
1541 " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n"
1542 " N = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V)) = "
1543 "%d (1 + 1 / %f (%d / %d - 1 - %d / %d)) = %f (%f)",
1544 Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta),
1545 Int_t(nAccepted), trigEff, Int_t(nTriggered),
1546 Int_t(nWithVertex), Int_t(beta), Int_t(nWithVertex),
1549 if (scheme & kTriggerEfficiency && !(scheme & kZeroBin)) {
1552 AliInfo(Form("Correcting for trigger efficiency:\n"
1553 " N = 1 / E_X * N = 1 / %f * %d = %f (%f)",
1554 trigEff, Int_t(ntotal), ntotal / trigEff, scaler));
1558 " Total of %9d events for %s\n"
1559 " of these %9d have an offline trigger\n"
1560 " of these N_T = %9d has the selected trigger\n"
1561 " of these N_V = %9d has a vertex\n"
1562 " of these N_A = %9d were in the selected range\n"
1563 " Triggers by hardware type:\n"
1565 " N_ac = %9d (%d+%d)\n"
1567 " Vertex efficiency: %f\n"
1568 " Trigger efficiency: %f\n"
1569 " Total number of events: N = %f\n"
1570 " Scaler (N_A/N): %f",
1571 Int_t(nAll), GetTitle(), Int_t(nOffline),
1572 Int_t(nTriggered), Int_t(nWithVertex), Int_t(nAccepted),
1573 Int_t(nB), Int_t(nA+nC), Int_t(nA), Int_t(nC), Int_t(nE),
1574 vtxEff, trigEff, ntotal, scaler));
1578 //________________________________________________________________________
1580 AliBasedNdetaTask::CentralityBin::GetResultName(Int_t rebin,
1582 const char* postfix) const
1585 n = Form("dndeta%s%s",GetName(), postfix);
1586 if (rebin > 1) n.Append(Form("_rebin%02d", rebin));
1587 if (sym) n.Append("_mirror");
1590 //________________________________________________________________________
1592 AliBasedNdetaTask::CentralityBin::GetResult(Int_t rebin,
1594 const char* postfix) const
1597 AliWarning(Form("No output list defined in %s [%3d,%3d]", GetName(),
1601 TString n = GetResultName(rebin, sym, postfix);
1602 TObject* o = fOutput->FindObject(n.Data());
1604 // AliWarning(Form("Object %s not found in output list", n.Data()));
1607 return static_cast<TH1*>(o);
1610 //________________________________________________________________________
1612 AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum,
1613 const char* postfix,
1616 const TH2F* shapeCorr,
1627 // Generate the dN/deta result from input
1630 // sum Sum of 2D hists
1631 // postfix Post fix on names
1632 // rootProj Whether to use ROOT TH2::ProjectionX
1633 // corrEmpty Correct for empty bins
1634 // shapeCorr Shape correction to use
1635 // scaler Event-level normalization scaler
1636 // symmetrice Whether to make symmetric extensions
1637 // rebin Whether to rebin
1638 // cutEdges Whether to cut edges when rebinning
1640 DGUARD(fDebug,1,"Make centrality bin result from %s", sum->GetName());
1641 TH2D* copy = static_cast<TH2D*>(sum->Clone(Form("d2Ndetadphi%s%s",
1642 GetName(), postfix)));
1643 TH1D* accNorm = ProjectX(sum, Form("norm%s%s",GetName(), postfix), 0, 0,
1644 rootProj, corrEmpty, false);
1645 accNorm->SetDirectory(0);
1647 // ---- Scale by shape correction ----------------------------------
1648 if (shapeCorr) copy->Divide(shapeCorr);
1649 else AliInfo("No shape correction specified, or disabled");
1651 // --- Normalize to the coverage -----------------------------------
1652 ScaleToCoverage(copy, accNorm);
1654 // --- Event-level normalization -----------------------------------
1655 copy->Scale(scaler);
1657 // --- Project on X axis -------------------------------------------
1658 TH1D* dndeta = ProjectX(copy, Form("dndeta%s%s",GetName(), postfix),
1659 1, sum->GetNbinsY(), rootProj, corrEmpty);
1660 dndeta->SetDirectory(0);
1661 // Event-level normalization
1662 dndeta->Scale(1., "width");
1663 copy->Scale(1., "width");
1665 TH1D* dndetaMCCorrection = 0;
1666 TList* centlist = 0;
1667 TH1D* dndetaMCtruth = 0;
1668 TList* truthcentlist = 0;
1670 // Possible final correction to <MC analysis> / <MC truth>
1672 centlist = static_cast<TList*> (mclist->FindObject(GetListName()));
1674 dndetaMCCorrection = static_cast<TH1D*>(centlist->FindObject(Form("dndeta%s%s",GetName(), postfix)));
1676 truthcentlist = static_cast<TList*> (truthlist->FindObject(GetListName()));
1678 dndetaMCtruth = static_cast<TH1D*> (truthcentlist->FindObject("dndetaTruth"));
1679 //std::cout<<dndetaMCCorrection<<" "<<dndetaMCtruth<<std::endl;
1680 if(dndetaMCCorrection && dndetaMCtruth) {
1681 AliInfo("Correcting with final MC correction");
1682 dndetaMCCorrection->Divide(dndetaMCtruth);
1683 dndeta->Divide(dndetaMCCorrection);
1685 //std::cout<<"histo "<<Form("dndeta%s%s",GetName(), postfix)<<" "<<GetListName()<<" "<<dndetaMCCorrection<<std::endl;
1686 //std::cout<<"truth "<<GetListName()<<" "<<dndetaMCtruth<<std::endl;
1689 else AliInfo("No final MC correction applied");
1691 // --- Set some histogram attributes -------------------------------
1693 Int_t rColor = GetColor(color);
1694 if (postfix && postfix[0] != '\0') post = Form(" (%s)", postfix);
1695 SetHistogramAttributes(dndeta, rColor, marker,
1696 Form("ALICE %s%s", GetName(), post.Data()));
1697 SetHistogramAttributes(accNorm, rColor, marker,
1698 Form("ALICE %s normalisation%s",
1699 GetName(), post.Data()));
1701 // --- Make symmetric extensions and rebinnings --------------------
1702 if (symmetrice) fOutput->Add(Symmetrice(dndeta));
1703 fOutput->Add(dndeta);
1704 fOutput->Add(accNorm);
1706 fOutput->Add(Rebin(dndeta, rebin, cutEdges));
1707 if (symmetrice) fOutput->Add(Symmetrice(Rebin(dndeta, rebin, cutEdges)));
1710 //________________________________________________________________________
1712 AliBasedNdetaTask::CentralityBin::End(TList* sums,
1715 const TH2F* shapeCorr,
1729 // End of processing
1732 // sums List of sums
1733 // results Output list of results
1734 // shapeCorr Shape correction or nil
1735 // trigEff Trigger efficiency
1736 // symmetrice Whether to symmetrice the results
1737 // rebin Whether to rebin the results
1738 // corrEmpty Whether to correct for empty bins
1739 // cutEdges Whether to cut edges when rebinning
1740 // triggerMask Trigger mask
1742 DGUARD(fDebug,1,"End centrality bin procesing");
1744 fSums = dynamic_cast<TList*>(sums->FindObject(GetListName()));
1746 AliError("Could not retrieve TList fSums");
1750 fOutput = new TList;
1751 fOutput->SetName(GetListName());
1752 fOutput->SetOwner();
1753 results->Add(fOutput);
1756 if (!ReadSum(fSums, false)) {
1757 AliInfo("This task did not produce any output");
1761 if (!fSumMC) ReadSum(fSums, true);
1763 fTriggers = static_cast<TH1I*>(fSums->FindObject("triggers"));
1766 AliError("Couldn't find histogram 'triggers' in list");
1770 // --- Get normalization scaler ------------------------------------
1771 Double_t epsilonT = trigEff;
1772 Double_t epsilonT0 = trigEff;
1774 if (triggerMask == AliAODForwardMult::kNSD) {
1775 // This is a local change
1777 AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT));
1779 else if (triggerMask == AliAODForwardMult::kInel) {
1780 // This is a local change
1782 AliWarning(Form("Using hard-coded Inel trigger efficiency of %f",epsilonT));
1784 if (scheme & kZeroBin) {
1785 if (triggerMask==AliAODForwardMult::kInel)
1786 epsilonT0 = 0.785021; // 0.100240;
1787 else if (triggerMask==AliAODForwardMult::kInelGt0)
1789 else if (triggerMask==AliAODForwardMult::kNSD)
1790 epsilonT0 = .706587;
1792 AliWarning(Form("Using hard-coded NCluster>0 trigger efficiency of %f",
1796 // Get our histograms
1798 TH2D* sum = fSum->CalcSum(fOutput, nSum, epsilonT0, 1,
1799 marker, rootProj, corrEmpty);
1800 Double_t nSumMC = 0;
1802 if (fSumMC) sumMC = fSumMC->CalcSum(fOutput, nSumMC,
1803 epsilonT0, 1, marker,
1804 rootProj, corrEmpty);
1806 AliError("Failed to get sum from summer - bailing out");
1810 Double_t ntotal = nSum;
1811 Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal);
1813 AliError("Failed to calculate normalization - bailing out");
1816 fOutput->Add(fTriggers->Clone());
1818 // --- Make result and store ---------------------------------------
1819 MakeResult(sum, "", rootProj, corrEmpty, (scheme & kShape) ? shapeCorr : 0,
1820 scaler, symmetrice, rebin, cutEdges, marker, color,
1823 // --- Process result from TrackRefs -------------------------------
1825 MakeResult(sumMC, "MC", rootProj, corrEmpty,
1826 (scheme & kShape) ? shapeCorr : 0,
1827 scaler, symmetrice, rebin, cutEdges,
1828 GetMarkerStyle(GetMarkerBits(marker)+4), color,
1832 // if (!IsAllBin()) return;