1 #include "AliForwardMultDists.h"
2 #include <AliForwardUtil.h>
3 #include <AliAODForwardMult.h>
4 #include <AliAODCentralMult.h>
5 #include <AliAODEvent.h>
11 #include <TParameter.h>
15 //____________________________________________________________________
16 AliForwardMultDists::AliForwardMultDists()
17 : AliAnalysisTaskSE(),
36 //____________________________________________________________________
37 AliForwardMultDists::AliForwardMultDists(const char* name)
38 : AliAnalysisTaskSE(name),
46 fTriggerMask(AliAODForwardMult::kV0AND),
59 * @param name Name of the task
61 DefineOutput(1, TList::Class());
62 DefineOutput(2, TList::Class());
65 //____________________________________________________________________
66 AliForwardMultDists::AliForwardMultDists(const AliForwardMultDists& o)
67 : AliAnalysisTaskSE(o),
69 fSymmetric(o.fSymmetric),
70 fNegative(o.fNegative),
71 fPositive(o.fPositive),
73 fTriggers(o.fTriggers),
75 fTriggerMask(o.fTriggerMask),
78 fFirstEvent(o.fFirstEvent),
79 fForwardCache(o.fForwardCache),
80 fCentralCache(o.fCentralCache),
83 fUsePhiAcc(o.fUsePhiAcc)
85 //____________________________________________________________________
87 AliForwardMultDists::operator=(const AliForwardMultDists& o)
89 if (&o == this) return *this;
91 fSymmetric = o.fSymmetric;
92 fNegative = o.fNegative;
93 fPositive = o.fPositive;
95 fTriggers = o.fTriggers;
97 fTriggerMask = o.fTriggerMask;
100 fFirstEvent = o.fFirstEvent;
101 fForwardCache = o.fForwardCache;
102 fCentralCache = o.fCentralCache;
103 fMCCache = o.fMCCache;
105 fUsePhiAcc = o.fUsePhiAcc;
110 //____________________________________________________________________
112 AliForwardMultDists::UserCreateOutputObjects()
115 * Create output objects - called at start of job in slave
120 fList->SetName("myMultAna");
123 fTriggers = AliAODForwardMult::MakeTriggerHistogram("triggers",
125 fTriggers->SetDirectory(0);
127 fStatus = AliAODForwardMult::MakeStatusHistogram("status");
128 fStatus->SetDirectory(0);
130 fList->Add(fTriggers);
135 //____________________________________________________________________
137 AliForwardMultDists::UserExec(Option_t* /*option=""*/)
140 * Analyse a single event
142 * @param option Not used
144 AliAODEvent* aod = AliForwardUtil::GetAODEvent(this);
147 TObject* obj = aod->FindListObject("Forward");
149 AliWarning("No forward object found");
152 AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
154 obj = aod->FindListObject("CentralClusters");
156 AliWarning("No central object found");
159 AliAODCentralMult* central = static_cast<AliAODCentralMult*>(obj);
162 obj = aod->FindListObject("primary");
163 // We should have a forward object at least
164 if (obj) primary = static_cast<TH2D*>(obj);
166 const TH2& forwardData = forward->GetHistogram();
167 const TH2& centralData = central->GetHistogram();
170 StoreInformation(forward);
171 SetupForData(forwardData, primary != 0);
175 if (!forward->CheckEvent(fTriggerMask, fMinIpZ, fMaxIpZ, 0, 0,
176 fTriggers, fStatus)) return;
178 // Info("", "Event vz=%f", forward->GetIpZ());
180 ProjectX(forwardData, *fForwardCache, fUsePhiAcc);
181 ProjectX(centralData, *fCentralCache, fUsePhiAcc);
182 ProjectX(primary, fMCCache);
186 while ((bin = static_cast<EtaBin*>(next()))) {
187 bin->Process(*fForwardCache, *fCentralCache,
188 forwardData, centralData,
204 kDownTriangle = 0x008,
210 * Get a ROOT marker style from bit options
212 * @param bits Bit mask of options
214 * @return ROOT marker style
216 Int_t GetMarkerStyle(UShort_t bits)
218 Int_t base = bits & (0xFE);
219 Bool_t hollow = bits & kHollow;
221 case kCircle: return (hollow ? 24 : 20);
222 case kSquare: return (hollow ? 25 : 21);
223 case kUpTriangle: return (hollow ? 26 : 22);
224 case kDownTriangle: return (hollow ? 32 : 23);
225 case kDiamond: return (hollow ? 27 : 33);
226 case kCross: return (hollow ? 28 : 34);
227 case kStar: return (hollow ? 30 : 29);
232 * Get the marker option bits from a ROOT style
234 * @param style ROOT style
236 * @return Pattern of marker options
238 UShort_t GetMarkerBits(Int_t style)
242 case 24: case 25: case 26: case 27: case 28: case 30: case 32:
243 bits |= kHollow; break;
246 case 20: case 24: bits |= kCircle; break;
247 case 21: case 25: bits |= kSquare; break;
248 case 22: case 26: bits |= kUpTriangle; break;
249 case 23: case 32: bits |= kDownTriangle; break;
250 case 27: case 33: bits |= kDiamond; break;
251 case 28: case 34: bits |= kCross; break;
252 case 29: case 30: bits |= kStar; break;
256 static Int_t GetIndexMarker(Int_t idx)
258 const UShort_t nMarkers = 7;
259 UShort_t markers[] = {
269 return markers[idx % nMarkers];
272 * Flip the 'hollow-ness' of a marker
274 * @param style ROOT Style
278 Int_t FlipHollowStyle(Int_t style)
280 UShort_t bits = GetMarkerBits(style);
281 Int_t ret = GetMarkerStyle(bits ^ kHollow);
285 //____________________________________________________________________
287 AliForwardMultDists::Terminate(Option_t* /*option=""*/)
290 * Called at the end of the final processing of the job on the
291 * full data set (merged data)
294 * @param option Not used
296 TList* out = new TList;
297 out->SetName("results");
300 TList* in = dynamic_cast<TList*>(GetOutputData(1));
302 AliError("No data connected to slot 1 here");
305 out->Add(in->FindObject("triggers")->Clone());
306 out->Add(in->FindObject("status")->Clone());
308 THStack* sym = new THStack("all", "All distributions");
309 THStack* pos = new THStack("all", "All distributions");
310 THStack* neg = new THStack("all", "All distributions");
311 THStack* oth = new THStack("all", "All distributions");
315 while ((bin = static_cast<EtaBin*>(next()))) {
316 bin->Terminate(in, out, fMaxN);
319 if (bin->IsSymmetric()) sta = sym;
320 else if (bin->IsNegative()) sta = neg;
321 else if (bin->IsPositive()) sta = pos;
323 TH1* hh[] = { bin->fSum, bin->fTruth, 0 };
326 Int_t idx = sta->GetHists() ? sta->GetHists()->GetEntries() : 0;
327 if (bin->fTruth) idx /= 2;
329 Int_t mkrBits = GetIndexMarker(idx);
331 Int_t factor = TMath::Power(10, idx);
332 TH1* h = static_cast<TH1*>((*ph)->Clone());
335 h->SetTitle(Form("%s (#times%d)", h->GetTitle(), Int_t(factor)));
336 h->SetMarkerStyle(GetMarkerStyle(mkrBits));
343 TList* lsym = static_cast<TList*>(out->FindObject("symmetric"));
344 TList* lneg = static_cast<TList*>(out->FindObject("negative"));
345 TList* lpos = static_cast<TList*>(out->FindObject("positive"));
346 TList* loth = static_cast<TList*>(out->FindObject("other"));
348 if (lsym) lsym->Add(sym);
349 if (lneg) lneg->Add(neg);
350 if (lpos) lpos->Add(pos);
351 if (loth) loth->Add(oth);
356 //____________________________________________________________________
358 AliForwardMultDists::StoreInformation(const AliAODForwardMult* aod)
360 fList->Add(AliForwardUtil::MakeParameter("sys", aod->GetSystem()));
361 fList->Add(AliForwardUtil::MakeParameter("snn", aod->GetSNN()));
362 fList->Add(AliForwardUtil::MakeParameter("trigger", ULong_t(fTriggerMask)));
363 fList->Add(AliForwardUtil::MakeParameter("minIpZ", fMinIpZ));
364 fList->Add(AliForwardUtil::MakeParameter("maxIpZ", fMaxIpZ));
365 fList->Add(AliForwardUtil::MakeParameter("maxN", UShort_t(fMaxN)));
368 //____________________________________________________________________
370 AliForwardMultDists::SetupForData(const TH2& hist, Bool_t useMC)
373 * Set-up internal structures on first seen event
375 * @param hist Basic histogram template from AOD object
377 TAxis* xaxis = hist.GetXaxis();
378 if (!xaxis->GetXbins() || xaxis->GetXbins()->GetSize() <= 0) {
379 fForwardCache = new TH1D("forwardCache", "Projection of Forward data",
380 xaxis->GetNbins(), xaxis->GetXmin(),
382 fCentralCache = new TH1D("centralCache", "Projection of Central data",
383 xaxis->GetNbins(), xaxis->GetXmin(),
386 fMCCache = new TH1D("mcCache", "Projection of Mc data",
387 xaxis->GetNbins(), xaxis->GetXmin(),
391 fForwardCache = new TH1D("forwardCache", "Projection of Forward data",
392 xaxis->GetNbins(),xaxis->GetXbins()->GetArray());
393 fCentralCache = new TH1D("centralCache", "Projection of Central data",
394 xaxis->GetNbins(),xaxis->GetXbins()->GetArray());
396 fMCCache = new TH1D("mcCache", "Projection of Mc data",
397 xaxis->GetNbins(),xaxis->GetXbins()->GetArray());
399 fForwardCache->SetDirectory(0);
400 fForwardCache->GetXaxis()->SetTitle("#eta");
401 fForwardCache->GetYaxis()->SetTitle("#sum#frac{d^{2}N}{d#etadphi}");
402 fForwardCache->Sumw2();
403 fCentralCache->SetDirectory(0);
404 fCentralCache->GetXaxis()->SetTitle("#eta");
405 fCentralCache->GetYaxis()->SetTitle("#sum#frac{d^{2}N}{d#etadphi}");
406 fCentralCache->Sumw2();
408 fMCCache->SetDirectory(0);
409 fMCCache->GetXaxis()->SetTitle("#eta");
410 fMCCache->GetYaxis()->SetTitle("#sum#frac{d^{2}N}{d#etadphi}");
416 while ((bin = static_cast<EtaBin*>(next()))) {
417 bin->SetupForData(fList, hist, fMaxN, useMC);
421 //____________________________________________________________________
423 AliForwardMultDists::ProjectX(const TH2& input, TH1& cache, Bool_t usePhiAcc)
426 * Project a 2D histogram into a 1D histogram taking care to use
427 * either the @f$\phi2f$ acceptance stored in the overflow bins, or
428 * the @f$\eta@f$ coverage stored in the underflow bins.
430 * @param input 2D histogram to project
431 * @param cache 1D histogram to project into
432 * @param usePhiAcc If true, use the @f$\phi2f$ acceptance stored in
433 * the overflow bins, or if false the @f$\eta@f$ coverage stored in
434 * the underflow bins.
438 Int_t nPhi = input.GetNbinsY();
439 Int_t nEta = input.GetNbinsX();
441 for (Int_t iEta = 1; iEta <= nEta; iEta++) {
442 Double_t phiAcc = input.GetBinContent(iEta, nPhi+1);
443 Double_t etaCov = input.GetBinContent(iEta, 0);
444 Double_t factor = usePhiAcc ? phiAcc : etaCov;
446 if (factor < 1e-3) continue;
449 for (Int_t iPhi = 1; iPhi <= nPhi; iPhi++) {
450 Double_t c = input.GetBinContent(iEta, iPhi);
451 Double_t e = input.GetBinError(iEta, iPhi);
456 cache.SetBinContent(iEta, factor * sum);
457 cache.SetBinError(iEta, factor * TMath::Sqrt(e2sum));
460 //____________________________________________________________________
462 AliForwardMultDists::ProjectX(const TH2* input, TH1* cache)
465 * Project on @f$\eta@f$ axis. If any of the pointers passed is
471 if (!input || !cache) return;
474 Int_t nPhi = input->GetNbinsY();
475 Int_t nEta = input->GetNbinsX();
477 for (Int_t iEta = 1; iEta <= nEta; iEta++) {
480 for (Int_t iPhi = 1; iPhi <= nPhi; iPhi++) {
481 Double_t c = input->GetBinContent(iEta, iPhi);
482 Double_t e = input->GetBinError(iEta, iPhi);
487 cache->SetBinContent(iEta, sum);
488 cache->SetBinError(iEta, TMath::Sqrt(e2sum));
491 //____________________________________________________________________
493 AliForwardMultDists::AddBin(Double_t etaLow, Double_t etaMax)
496 * Add an @f$\eta@f$ bin
498 * @param etaLow Low cut on @f$\eta@f$
499 * @param etaMax High cut on @f$\eta@f$
502 if (etaMax >= kInvalidEta) {
503 etaLow = -TMath::Abs(etaLow);
504 etaMax = +TMath::Abs(etaLow);
506 EtaBin* bin = new EtaBin(etaLow, etaMax);
507 // AliInfoF("Adding bin %f, %f: %s", etaLow, etaMax, bin->GetName());
510 //____________________________________________________________________
512 AliForwardMultDists::SetTriggerMask(const char* mask)
515 * Set the trigger mask
519 fTriggerMask = AliAODForwardMult::MakeTriggerMask(mask);
521 //____________________________________________________________________
523 AliForwardMultDists::Print(Option_t* /*option=""*/) const
528 * @param option Not used
530 std::cout << "Task: " << GetName() << " " << ClassName() << "\n"
532 << AliAODForwardMult::GetTriggerString(fTriggerMask) << "\n"
533 << " IpZ range: " << fMinIpZ <<" - "<< fMaxIpZ << "\n"
534 << " Max Nch: " << fMaxN << "\n"
535 << " Use phi acceptance: " << fUsePhiAcc << "\n"
536 << " Bins:" << std::endl;
539 while ((bin = static_cast<EtaBin*>(next()))) {
540 std::cout << " " << bin->GetName() << std::endl;
544 //====================================================================
545 AliForwardMultDists::EtaBin::EtaBin()
562 //____________________________________________________________________
563 AliForwardMultDists::EtaBin::EtaBin(Double_t minEta, Double_t maxEta)
578 * @param minEta Least @f$\eta@f$ to consider
579 * @param maxEta Largest @f$\eta@f$ to consider
581 fName = TString::Format("%+05.2f_%+05.2f", fMinEta, fMaxEta);
582 fName.ReplaceAll("-", "m");
583 fName.ReplaceAll("+", "p");
584 fName.ReplaceAll(".", "d");
586 //____________________________________________________________________
587 AliForwardMultDists::EtaBin::EtaBin(const EtaBin& o)
596 fResponse(o.fResponse),
598 fCoverage(o.fCoverage)
600 //____________________________________________________________________
601 AliForwardMultDists::EtaBin&
602 AliForwardMultDists::EtaBin::operator=(const EtaBin& o)
604 if (&o == this) return *this;
613 fResponse = o.fResponse;
615 fCoverage = o.fCoverage;
619 //____________________________________________________________________
621 AliForwardMultDists::EtaBin::IsSymmetric() const
623 return TMath::Abs(fMaxEta + fMinEta) < 1e-6;
625 //____________________________________________________________________
627 AliForwardMultDists::EtaBin::IsNegative() const
629 return TMath::Abs(fMaxEta) < 1e-6 && fMinEta < 0;
631 //____________________________________________________________________
633 AliForwardMultDists::EtaBin::IsPositive() const
635 return TMath::Abs(fMinEta) < 1e-6 && fMaxEta > 0;
637 //____________________________________________________________________
639 AliForwardMultDists::EtaBin::ParentName() const
641 if (IsSymmetric()) return "symmetric";
642 else if (IsNegative()) return "negative";
643 else if (IsPositive()) return "positive";
646 //____________________________________________________________________
648 AliForwardMultDists::EtaBin::FindParent(TList* l, Bool_t create) const
650 const char* parent = ParentName();
651 TObject* op = l->FindObject(parent);
653 if (op) return static_cast<TList*>(op);
654 if (!create) return 0;
656 // Info("FindParent", "Parent %s not found in %s, creating and adding",
657 // parent, l->GetName());
658 TList* p = new TList;
663 TList* lowEdges = new TList;
664 lowEdges->SetName("lowEdges");
665 lowEdges->SetOwner();
668 TList* highEdges = new TList;
669 highEdges->SetName("highEdges");
670 highEdges->SetOwner();
676 //____________________________________________________________________
678 AliForwardMultDists::EtaBin::SetupForData(TList* list, const TH2& hist,
679 UShort_t max, Bool_t useMC)
682 * Set-up internal structures on first event.
684 * @param list List to add information to
685 * @param hist Template histogram
686 * @param max Maximum number of particles
687 * @param useMC Whether to set-up for MC input
689 TList* p = FindParent(list, true);
690 TList* l = new TList;
691 l->SetName(GetName());
694 TList* le = static_cast<TList*>(p->FindObject("lowEdges"));
695 TList* he = static_cast<TList*>(p->FindObject("highEdges"));
697 AliError("Failed to get bin edge lists from parent");
701 Int_t n = le->GetEntries();
702 TParameter<double>* lp =
703 new TParameter<double>(Form("minEta%02d", n), fMinEta);
704 TParameter<double>* hp =
705 new TParameter<double>(Form("maxEta%02d", n), fMaxEta);
706 lp->SetMergeMode('f');
707 hp->SetMergeMode('f');
712 fMinBin = hist.GetXaxis()->FindBin(fMinEta);
713 fMaxBin = hist.GetXaxis()->FindBin(fMaxEta-.00001);
715 fSum = new TH1D("rawDist",
716 Form("Raw P(#it{N}_{ch}) in %+5.2f<#eta<%+5.2f",
719 fSum->SetDirectory(0);
720 fSum->GetXaxis()->SetTitle("#it{N}_{ch}");
721 fSum->GetYaxis()->SetTitle("Raw P(#it{N}_{ch})");
722 fSum->SetFillColor(kRed+1);
723 fSum->SetFillStyle(0);
724 // fSum->SetOption("hist e");
725 fSum->SetMarkerStyle(20);
726 fSum->SetMarkerColor(kRed+1);
727 fSum->SetLineColor(kBlack);
731 fCorr = new TH2D("corr", "Correlation of SPD and FMD signals",
732 max+2, -1.5, max+.5, max+2, -1.5, max+.5);
733 fCorr->SetDirectory(0);
734 fCorr->GetXaxis()->SetTitle("Forward");
735 fCorr->GetYaxis()->SetTitle("Central");
736 fCorr->SetOption("colz");
739 fCoverage = new TH1D("coverage", "Fraction of bins with coverage",
741 fCoverage->SetDirectory(0);
742 fCoverage->SetXTitle("Fraction of bins [%]");
743 fCoverage->SetFillStyle(3001);
744 fCoverage->SetFillColor(kGreen+1);
748 fResponse = new TH2D("response", "Reponse matrix",
749 max+1, -.5, max+.5, max+1, -.5, max+.5);
750 fResponse->SetDirectory(0);
751 fResponse->SetXTitle("MC");
752 fResponse->SetYTitle("Signal");
753 fResponse->SetOption("colz");
756 fTruth = static_cast<TH1D*>(fSum->Clone("truth"));
757 fTruth->SetTitle(Form("True P(#it{N}_{ch}) in %+5.2f<#eta<%+5.2f",
759 fTruth->SetLineColor(kBlack);
760 fTruth->SetFillColor(kBlue+1);
761 fTruth->SetFillStyle(0);
762 fTruth->SetDirectory(0);
763 /// fTruth->SetOption("");
764 fTruth->SetMarkerColor(kBlue+1);
765 fTruth->SetMarkerStyle(24);
769 //____________________________________________________________________
771 AliForwardMultDists::EtaBin::Process(const TH1& sumForward,
772 const TH1& sumCentral,
778 * Process a single event
780 * @param sumForward Projection of forward data
781 * @param sumCentral Projection of the central data
782 * @param forward The original forward data
783 * @param central The original central data
791 Double_t mce2sum = 0;
792 for (Int_t iEta = fMinBin; iEta <= fMaxBin; iEta++) {
793 Double_t cF = sumForward.GetBinContent(iEta);
794 Double_t eF = sumForward.GetBinError(iEta);
795 Bool_t bF = forward.GetBinContent(iEta, 0) > 0;
796 Double_t cC = sumCentral.GetBinContent(iEta);
797 Double_t eC = sumCentral.GetBinError(iEta);
798 Bool_t bC = central.GetBinContent(iEta, 0) > 0;
799 Double_t cM = (mc ? mc->GetBinContent(iEta) : 0);
800 Double_t eM = (mc ? mc->GetBinError(iEta) : 0);
804 "bF=%d cF=%7.4f+/-%8.5f "
805 "bC=%d cC=%7.4f+/-%8.5f "
807 bF, cF, eF, bC, cC, eC, cM, eM);
816 // If we have an overlap - as given by the eta-coverage,
817 // calculate the mean
820 e = TMath::Sqrt(eF*eF + eC*eC);
823 // Else, if we have coverage from forward, use that
829 // Else, if we have covrage from central, use that
835 // Else, we have incomplete coverage
838 if (fsum < 0) fsum = 0;
842 if (csum < 0) csum = 0;
851 Double_t w = 1; // e2sum > 0 ? 1/TMath::Sqrt(e2sum) : 1
853 fCorr->Fill(fsum, csum);
855 Int_t nTotal = fMaxBin - fMinBin + 1;
856 fCoverage->Fill(100*float(covered) / nTotal);
857 // Info(GetName(), "covered: %3d, nTotal: %3d -> %3d%% mc: %3d",
858 // covered, nTotal, int(100*float(covered)/nTotal), int(mcsum));
862 w = 1; // mce2sum > 0 ? 1/TMath::Sqrt(mce2sum) : 1
863 fTruth->Fill(mcsum, w);
866 fResponse->Fill(mcsum, sum);
869 //____________________________________________________________________
871 AliForwardMultDists::EtaBin::Terminate(TList* in, TList* out, UShort_t maxN)
874 * Called at the end of the final processing of the job on the
875 * full data set (merged data)
877 * @param in Input list
878 * @param out Output list
879 * @param maxN Maximum number of @f$N_{ch}@f$ to consider
881 TList* ip = FindParent(in, false);
883 AliErrorF("Parent folder %s not found in input", ParentName());
887 TList* i = dynamic_cast<TList*>(ip->FindObject(GetName()));
889 AliErrorF("List %s not found in input", GetName());
893 TList* op = FindParent(out, true);
894 TList* o = static_cast<TList*>(i->Clone());
898 fSum = static_cast<TH1*>(o->FindObject("rawDist"));
899 fTruth = static_cast<TH1*>(o->FindObject("truth"));
900 TH1* hists[] = { fSum, fTruth, 0 };
905 Double_t intg = h->Integral(1, maxN);
911 //====================================================================