TList* mclist = 0;
TList* truthlist = 0;
- if(fFinalMCCorrFile.Contains(".root")) {
+ if (fFinalMCCorrFile.Contains(".root")) {
TFile* ftest = TFile::Open(fFinalMCCorrFile.Data());
if(ftest) {
- mclist = dynamic_cast<TList*> (ftest->Get(Form("%sResults", GetName())));
- truthlist = dynamic_cast<TList*> (ftest->Get("MCTruthResults"));
+ mclist = dynamic_cast<TList*>(ftest->Get(Form("%sResults",GetName())));
+ truthlist = dynamic_cast<TList*>(ftest->Get("MCTruthResults"));
}
- else AliWarning("MC analysis file invalid - no final MC correction possible");
-
+ else
+ AliWarning("MC analysis file invalid - no final MC correction possible");
}
Int_t style = GetMarker();
Int_t color = GetColor();
if(mclist)
centlist = static_cast<TList*> (mclist->FindObject(GetListName()));
if(centlist)
- dndetaMCCorrection = static_cast<TH1D*> (centlist->FindObject(Form("dndeta%s%s",GetName(), postfix)));
+ dndetaMCCorrection = static_cast<TH1D*>(centlist->FindObject(Form("dndeta%s%s",GetName(), postfix)));
if(truthlist)
truthcentlist = static_cast<TList*> (truthlist->FindObject(GetListName()));
if(truthcentlist)
// --- Make result and store ---------------------------------------
MakeResult(sum, "", rootProj, corrEmpty, (scheme & kShape) ? shapeCorr : 0,
- scaler, symmetrice, rebin, cutEdges, marker, color, mclist, truthlist);
+ scaler, symmetrice, rebin, cutEdges, marker, color,
+ mclist, truthlist);
// --- Process result from TrackRefs -------------------------------
if (sumMC)
MakeResult(sumMC, "MC", rootProj, corrEmpty,
(scheme & kShape) ? shapeCorr : 0,
scaler, symmetrice, rebin, cutEdges,
- GetMarkerStyle(GetMarkerBits(marker)+4), color, mclist, truthlist);
+ GetMarkerStyle(GetMarkerBits(marker)+4), color,
+ mclist, truthlist);
// Temporary stuff
// if (!IsAllBin()) return;
#include <TLatex.h>
#include <TImage.h>
+Double_t myFunc(Double_t* xp, Double_t* pp);
+
/**
* Class to draw dN/deta results
*
*
*/
dNdetaDrawer()
- : fShowOthers(false), // Bool_t
+ : fShowOthers(0), // Bool_t
fShowAlice(false), // Bool_t
fShowRatios(false), // Bool_t
fShowLeftRight(false), // Bool_t
fCentAxis(0), // TAxis*
fTriggers(0), // TH1*
fTruth(0), // TH1*
- fRangeParam(0)
-
+ fRangeParam(0),
+ fFwdSysErr(0.076),
+ fCenSysErr(0),
+ fRemoveOuters(false),
+ fClusterScale("")
{
fRangeParam = new RangeParam;
fRangeParam->fMasterAxis = 0;
*
* @param x Whether to show or not
*/
- void SetShowOthers(Bool_t x) { fShowOthers = x; }
+ void SetShowOthers(UShort_t x) { fShowOthers = x; }
//__________________________________________________________________
/**
* Show ALICE published data
* @param x To show or not
*/
void SetShowLeftRight(Bool_t x) { fShowLeftRight = x; }
+ //__________________________________________________________________
+ /**
+ * Whether to show rings
+ *
+ * @param x To show or not
+ */
void SetShowRings(Bool_t x) { fShowRings = x; }
//__________________________________________________________________
/**
* @param x Title
*/
void SetTitle(TString x) { fTitle = x; }
+ //__________________________________________________________________
+ /**
+ * Set the systematic error in the forward region
+ *
+ * @param e Systematic error in the forward region
+ */
+ void SetForwardSysError(Double_t e=0) { fFwdSysErr = e; }
+ //__________________________________________________________________
+ /**
+ * Set the systematic error in the forward region
+ *
+ * @param e Systematic error in the forward region
+ */
+ void SetCentralSysError(Double_t e=0) { fCenSysErr = e; }
/* @} */
//==================================================================
/**
fOthers = new TMultiGraph();
// --- Loop over input data --------------------------------------
- FetchResults(mcTruth, "MCTruth", max, rmax, amax);
- FetchResults(forward, "Forward", max, rmax, amax);
- FetchResults(clusters, "Central", max, rmax, amax);
+ /*TObjArray* mcA =*/ FetchResults(mcTruth, "MCTruth", max, rmax, amax);
+ TObjArray* fwdA = FetchResults(forward, "Forward", max, rmax, amax);
+ TObjArray* cenA = FetchResults(clusters, "Central", max, rmax, amax);
// --- Get trigger information -----------------------------------
TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
if (!fOthers->GetListOfGraphs() ||
fOthers->GetListOfGraphs()->GetEntries() <= 0) {
Warning("Run", "No other data found - disabling that");
- fShowOthers = false;
+ fShowOthers = 0;
}
if (!fRatios->GetHists() ||
fRatios->GetHists()->GetEntries() <= 0) {
// fLeftRight->ls();
fShowLeftRight = false;
}
+ if (fFwdSysErr > 0) {
+ if (fCenSysErr <= 0) fCenSysErr = fFwdSysErr;
+ for (Int_t i = 0; i < fwdA->GetEntriesFast(); i++) {
+ TH1* fwd = static_cast<TH1*>(fwdA->At(i));
+ TH1* cen = static_cast<TH1*>(cenA->At(i));
+ CorrectForward(fwd);
+ CorrectCentral(cen);
+ Double_t low, high;
+ TH1* tmp = Merge(cen, fwd, low, high);
+ TF1* f = FitMerged(tmp, low, high);
+ MakeSysError(tmp, cen, fwd, f);
+ delete f;
+ Info("", "Adding systematic error histogram %s",
+ tmp->GetName());
+ fResults->GetHists()->AddFirst(tmp, "e5");
+ TH1* tmp2 = Symmetrice(tmp);
+ tmp2->SetFillColor(tmp->GetFillColor());
+ tmp2->SetFillStyle(tmp->GetFillStyle());
+ tmp2->SetMarkerStyle(tmp->GetMarkerStyle());
+ tmp2->SetLineWidth(tmp->GetLineWidth());
+ fResults->GetHists()->AddFirst(tmp2, "e5");
+ fResults->Modified();
+ }
+ }
+ delete fwdA;
+ delete cenA;
// --- Close the input file --------------------------------------
file->Close();
TMultiGraph* FetchOthers(UShort_t centLow, UShort_t centHigh)
{
TMultiGraph* thisOther = 0;
- if (!fShowOthers && !fShowAlice) return 0;
+ if (fShowOthers == 0) return 0;
- Bool_t onlya = (fShowOthers ? false : true);
- Bool_t nomc = true;
UShort_t sys = (fSysString ? fSysString->GetUniqueID() : 0);
UShort_t trg = (fTrigString ? fTrigString->GetUniqueID() : 0);
UShort_t snn = (fSNNString ? fSNNString->GetUniqueID() : 0);
- Long_t ret = gROOT->ProcessLine(Form("GetData(%d,%d,%d,%d,%d,%d,%d);",
+ Long_t ret = gROOT->ProcessLine(Form("GetData(%d,%d,%d,%d,%d,%d);",
sys,snn,trg,
- centLow,centHigh,onlya,nomc));
- if (!ret) {
-#if 0
- Warning("FetchOthers", "No other data found for sys=%d, sNN=%d, "
- "trigger=%d %d%%-%d%% central %s",
- sys, snn, trg, centLow, centHigh,
- onlya ? " (ALICE results)" : "all");
-#endif
- return 0;
- }
- thisOther = reinterpret_cast<TMultiGraph*>(ret);
-
+ centLow,centHigh,
+ fShowOthers));
+ if (!ret) return 0;
+
+ thisOther = reinterpret_cast<TMultiGraph*>(ret);
return thisOther;
}
//__________________________________________________________________
* @param rmax On return, maximum of ratios
* @param amax On return, maximum of left-right comparisons
*/
- void FetchResults(const TList* list,
- const char* name,
- Double_t& max,
- Double_t& rmax,
- Double_t& amax)
+ TObjArray*
+ FetchResults(const TList* list,
+ const char* name,
+ Double_t& max,
+ Double_t& rmax,
+ Double_t& amax)
{
- UShort_t n = fCentAxis ? fCentAxis->GetNbins() : 0;
+ UShort_t n = fCentAxis ? fCentAxis->GetNbins() : 0;
if (n == 0) {
TList* all = static_cast<TList*>(list->FindObject("all"));
- if (!all)
+ if (!all) {
Error("FetchResults", "Couldn't find list 'all' in %s",
list->GetName());
- else
- FetchResults(all, name, FetchOthers(0,0), -1000, 0, max, rmax, amax);
- return;
+ return 0;
+ }
+ TObjArray* a = new TObjArray;
+ TH1* h = FetchResults(all, name, FetchOthers(0,0),
+ -1000, 0, max, rmax, amax);
+ a->AddAt(h, 0);
+ return a;
}
+ TObjArray* a = new TObjArray;
for (UShort_t i = 0; i < n; i++) {
UShort_t centLow = fCentAxis->GetBinLowEdge(i+1);
UShort_t centHigh = fCentAxis->GetBinUpEdge(i+1);
Int_t col = GetCentralityColor(i+1);
TString centTxt = Form("%3d%%-%3d%% central", centLow, centHigh);
- if (!thisCent)
+ if (!thisCent) {
Error("FetchResults", "Couldn't find list '%s' in %s",
lname.Data(), list->GetName());
- else
- FetchResults(thisCent, name, FetchOthers(centLow, centHigh),
- col, centTxt.Data(), max, rmax, amax);
+ continue;
+ }
+ TH1* h = FetchResults(thisCent, name, FetchOthers(centLow, centHigh),
+ col, centTxt.Data(), max, rmax, amax);
+ a->AddAt(h, i);
}
+ return a;
}
//__________________________________________________________________
Int_t GetCentralityColor(Int_t bin) const
* @param rmax On return, ratio maximum
* @param amax On return, left-right maximum
*/
- void FetchResults(const TList* list,
+ TH1* FetchResults(const TList* list,
const char* name,
TMultiGraph* thisOther,
Int_t color,
TH1* dndetaTruth = FetchResult(list, "dndetaTruth");
TH1* dndetaSym = 0;
TH1* dndetaMCSym = 0;
- TH1* tracks = FetchResult(list, "tracks");
- if (tracks) tracks->SetTitle("ALICE Tracks");
SetAttributes(dndeta, color);
SetAttributes(dndetaMC, fCentAxis ? color : color+2);
SetAttributes(dndetaTruth,color);
SetAttributes(dndetaSym, color);
SetAttributes(dndetaMCSym,fCentAxis ? color : color+2);
- SetAttributes(tracks, fCentAxis ? color : color+3);
if (dndetaMC && fCentAxis)
dndetaMC->SetMarkerStyle(dndetaMC->GetMarkerStyle()+2);
if (dndetaMCSym && fCentAxis)
ModifyTitle(dndetaTruth,centTxt);
ModifyTitle(dndetaSym, centTxt);
ModifyTitle(dndetaMCSym,centTxt);
- ModifyTitle(tracks, centTxt);
max = TMath::Max(max, AddHistogram(fResults, dndetaTruth, "e5 p"));
max = TMath::Max(max, AddHistogram(fResults, dndetaMC, dndetaMCSym));
max = TMath::Max(max, AddHistogram(fResults, dndeta, dndetaSym));
- max = TMath::Max(max, AddHistogram(fResults, tracks));
if (dndetaTruth) {
fTruth = dndetaTruth;
if (fShowLeftRight) {
fLeftRight->Add(Asymmetry(dndeta, amax));
fLeftRight->Add(Asymmetry(dndetaMC, amax));
- fLeftRight->Add(Asymmetry(tracks, amax));
}
if (thisOther) {
}
// fOthers->Add(thisOther);
}
- if (tracks) {
- if (!fRatios->GetHists() ||
- !fRatios->GetHists()->FindObject(tracks->GetName()))
- fRatios->Add(Ratio(dndeta, tracks, rmax));
- }
}
if (dndetaMC) {
fRatios->Add(Ratio(dndeta, dndetaMC, rmax));
fRatios->Add(Ratio(dndeta, fTruth, rmax));
fRatios->Add(Ratio(dndetaSym, fTruth, rmax));
}
+ return dndeta;
}
//__________________________________________________________________
/**
Double_t y1 = 0;
Double_t y2 = 0;
Double_t y3 = 0;
- if (!fShowRatios) w *= 1.4;
+ if (!fShowRatios) w *= 1.3;
else y1 = 0.3;
- if (!fShowLeftRight) w *= 1.4;
+ if (!fShowLeftRight) w *= 1.3;
else {
Double_t y11 = y1;
y1 = (y11 > 0.0001 ? 0.4 : 0.2);
l->SetBorderSize(0);
l->SetTextFont(132);
+ // Loop over items in stack and get unique items, while ignoring
+ // mirrored data and systematic error bands
TIter next(stack->GetHists());
TH1* hist = 0;
TObjArray unique;
unique.SetOwner();
+ Bool_t sysErrSeen = false;
while ((hist = static_cast<TH1*>(next()))) {
- TString n(hist->GetTitle());
- if (n.Contains("mirrored")) continue;
- if (unique.FindObject(n.Data())) continue;
- TObjString* s = new TObjString(n);
+ TString t(hist->GetTitle());
+ TString n(hist->GetName());
+ n.ToLower();
+ if (t.Contains("mirrored")) continue;
+ if (n.Contains("syserror")) { sysErrSeen = true; continue; }
+ if (unique.FindObject(t.Data())) continue;
+ TObjString* s = new TObjString(hist->GetTitle());
s->SetUniqueID(((hist->GetMarkerStyle() & 0xFFFF) << 16) |
((hist->GetMarkerColor() & 0xFFFF) << 0));
unique.Add(s);
// l->AddEntry(hist, hist->GetTitle(), "pl");
}
if (mg) {
+ // If we have other stuff, scan for unique names
TIter nexto(mg->GetListOfGraphs());
TGraph* g = 0;
while ((g = static_cast<TGraph*>(nexto()))) {
// l->AddEntry(hist, hist->GetTitle(), "pl");
}
}
- // unique.ls();
+
+ // Add legend entries for unique items only
TIter nextu(&unique);
TObject* s = 0;
Int_t i = 0;
else dd->SetMarkerColor(color);
dd->SetMarkerStyle(style);
}
+ if (sysErrSeen) {
+ // Add entry for systematic errors
+ TLegendEntry* d0 = l->AddEntry("d0", Form("%4.1f%% Systematic error",
+ 100*fFwdSysErr), "f");
+ d0->SetLineColor(kGray);
+ d0->SetMarkerColor(kGray);
+ d0->SetFillColor(kGray);
+ d0->SetFillStyle(3001);
+ d0->SetMarkerStyle(0);
+ d0->SetLineWidth(0);
+ i++;
+ }
+ if (!fCentAxis && i % 2 == 1) {
+ // To make sure the 'data' and 'mirrored' entries are on a line
+ // by themselves
+ TLegendEntry* dd = l->AddEntry("dd", " ", "");
+ dd->SetTextSize(0);
+ dd->SetFillStyle(0);
+ dd->SetFillColor(0);
+ dd->SetLineWidth(0);
+ dd->SetLineColor(0);
+ dd->SetMarkerSize(0);
+ }
+ // Add entry for 'data'
TLegendEntry* d1 = l->AddEntry("d1", "Data", "lp");
d1->SetLineColor(kBlack);
d1->SetMarkerColor(kBlack);
d1->SetMarkerStyle(20);
+ // Add entry for 'mirrored data'
TLegendEntry* d2 = l->AddEntry("d2", "Mirrored data", "lp");
d2->SetLineColor(kBlack);
d2->SetMarkerColor(kBlack);
p1->SetBorderSize(0);
p1->SetBorderMode(0);
p1->SetBottomMargin(yd > 0.001 ? 0.001 : 0.1);
- p1->SetRightMargin(0.05);
- p1->SetGridx();
+ p1->SetRightMargin(0.03);
+ if (fShowLeftRight || fShowRatios) p1->SetGridx();
p1->SetTicks(1,1);
p1->SetNumber(1);
p1->Draw();
fResults->SetMaximum(1.15*max);
fResults->SetMinimum(yd > 0.00001 ? -0.1 : 0);
- Double_t s = (yd > 0.00001 ? 1/(1-yd)/1.7 : 1.2);
- FixAxis(fResults, s, "#frac{1}{N} #frac{dN_{ch}}{d#eta}");
+ FixAxis(fResults, (1-yd)*(yd > .001 ? 1 : .9 / 1.2),
+ "#font[12]{#frac{1}{N} "
+ "#frac{dN_{#font[132]{ch}}}{d#font[152]{#eta}}}");
p1->Clear();
fResults->DrawClone("nostack e1");
fRangeParam->fSlave1Pad = p1;
// Draw other data
- if (fShowOthers || fShowAlice) {
+ if (fShowOthers != 0) {
TGraphAsymmErrors* o = 0;
TIter nextG(fOthers->GetListOfGraphs());
while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
logo++;
continue;
}
- TPad* pad = new TPad("logo", "logo", .12, .65, .25, .85, 0, 0, 0);
+ TPad* pad = new TPad("logo", "logo", .12, .7, .25, .9, 0, 0, 0);
pad->SetFillStyle(0);
pad->Draw();
pad->cd();
*/
void PlotRatios(Double_t max, Double_t y1, Double_t y2)
{
- if (!fShowRatios) return;
+ if (fShowRatios == 0) return;
bool isBottom = (y1 < 0.0001);
Double_t yd = y2 - y1;
// Make a sub-pad for the result itself
TPad* p2 = new TPad("p2", "p2", 0, y1, 1.0, y2, 0, 0, 0);
p2->SetTopMargin(0.001);
- p2->SetRightMargin(0.05);
+ p2->SetRightMargin(0.03);
p2->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
p2->SetGridx();
p2->SetTicks(1,1);
p2->cd();
// Fix up axis
- FixAxis(fRatios, 1/yd/1.7, "Ratios", 7);
+ FixAxis(fRatios, yd, "Ratios", 7);
fRatios->SetMaximum(1+TMath::Max(.22,1.05*max));
fRatios->SetMinimum(1-TMath::Max(.32,1.05*max));
// Make a sub-pad for the result itself
TPad* p3 = new TPad("p3", "p3", 0, y1, 1.0, y2, 0, 0, 0);
p3->SetTopMargin(0.001);
- p3->SetRightMargin(0.05);
+ p3->SetRightMargin(0.03);
p3->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
p3->SetGridx();
p3->SetTicks(1,1);
p3->cd();
// Fix up axis
- FixAxis(fLeftRight, 1/yd/1.7, "Right/Left", 4);
+ FixAxis(fLeftRight, yd, "Right/Left", 4);
fLeftRight->SetMaximum(1+TMath::Max(.12,1.05*max));
fLeftRight->SetMinimum(1-TMath::Max(.15,1.05*max));
* @param force Whether to draw the stack first or not
* @param ynDiv Divisions on Y axis
*/
- void FixAxis(THStack* stack, Double_t s, const char* ytitle,
+ void FixAxis(THStack* stack, Double_t yd, const char* ytitle,
Int_t ynDiv=210, Bool_t force=true)
{
if (!stack) {
Warning("FixAxis", "Stack %s has no histogram", stack->GetName());
return;
}
-
- h->SetXTitle("#eta");
+ Double_t s = 1/yd/1.2;
+ // Info("FixAxis", "for %s, s=1/%f=%f", stack->GetName(), yd, s);
+
+ h->SetXTitle("#font[152]{#eta}");
h->SetYTitle(ytitle);
TAxis* xa = h->GetXaxis();
TAxis* ya = h->GetYaxis();
+
+ // Int_t npixels = 20
+ // Float_t dy = gPad->AbsPixeltoY(0) - gPad->AbsPixeltoY(npixels);
+ // Float_t ts = dy/(gPad->GetY2() - gPad->GetY1());
+
if (xa) {
- xa->SetTitle("#eta");
+ // xa->SetTitle(h->GetXTitle());
// xa->SetTicks("+-");
xa->SetTitleSize(s*xa->GetTitleSize());
xa->SetLabelSize(s*xa->GetLabelSize());
xa->SetTickLength(s*xa->GetTickLength());
+ // xa->SetTitleOffset(xa->GetTitleOffset()/s);
if (stack != fResults) {
TAxis* rxa = fResults->GetXaxis();
}
}
if (ya) {
- ya->SetTitle(ytitle);
+ // ya->SetTitle(h->GetYTitle());
ya->SetDecimals();
// ya->SetTicks("+-");
ya->SetNdivisions(ynDiv);
ya->SetLabelSize(s*ya->GetLabelSize());
}
}
+ //__________________________________________________________________
+ /**
+ * Merge two histograms into one
+ *
+ * @param cen Central part
+ * @param fwd Forward part
+ * @param xlow On return, lower eta bound
+ * @param xhigh On return, upper eta bound
+ *
+ * @return Newly allocated histogram or null
+ */
+ TH1*
+ Merge(const TH1* cen, const TH1* fwd, Double_t& xlow, Double_t& xhigh)
+ {
+ TH1* tmp = static_cast<TH1*>(fwd->Clone("tmp"));
+ TString name(fwd->GetName());
+ name.ReplaceAll("Forward", "Merged");
+ tmp->SetName(name);
+
+ // tmp->SetMarkerStyle(28);
+ // tmp->SetMarkerColor(kBlack);
+ tmp->SetDirectory(0);
+ xlow = 100;
+ xhigh = -100;
+ for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
+ Double_t cc = cen->GetBinContent(i);
+ Double_t cf = fwd->GetBinContent(i);
+ Double_t ec = cen->GetBinError(i);
+ Double_t ef = fwd->GetBinError(i);
+ Double_t nc = cf;
+ Double_t ne = ef;
+ if (cc < 0.001 && cf < 0.01) continue;
+ xlow = TMath::Min(tmp->GetXaxis()->GetBinLowEdge(i),xlow);
+ xhigh = TMath::Max(tmp->GetXaxis()->GetBinUpEdge(i),xhigh);
+ if (cc > 0.001) {
+ nc = cc;
+ ne = ec;
+ if (cf > 0.001) {
+ nc = (cf + cc) / 2;
+ ne = TMath::Sqrt(ec*ec + ef*ef);
+ }
+ }
+ tmp->SetBinContent(i, nc);
+ tmp->SetBinError(i, ne);
+ }
+ return tmp;
+ }
+ //____________________________________________________________________
+ /**
+ * Fit @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$ to histogram data
+ *
+ * @param tmp Histogram
+ * @param xlow Lower x bound
+ * @param xhigh Upper x bound
+ *
+ * @return Fitted function
+ */
+ TF1*
+ FitMerged(TH1* tmp, Double_t xlow, Double_t xhigh)
+ {
+ TF1* tmpf = new TF1("tmpf", "gaus", xlow, xhigh);
+ tmp->Fit(tmpf, "NQ", "");
+ tmp->SetDirectory(0);
+
+ TF1* fit = new TF1("f", myFunc, xlow, xhigh, 4);
+ fit->SetParNames("a_{1}", "a_{2}", "#sigma_{1}", "#sigma_{2}");
+ fit->SetParameters(tmpf->GetParameter(0),
+ .2,
+ tmpf->GetParameter(2),
+ tmpf->GetParameter(2)/4);
+ fit->SetParLimits(3, 0, 100);
+ fit->SetParLimits(4, 0, 100);
+ tmp->Fit(fit,"0W","");
+
+ delete tmpf;
+ return fit;
+ }
+ //____________________________________________________________________
+ /**
+ * Make band of systematic errors
+ *
+ * @param tmp Histogram
+ * @param fit Fit
+ */
+ void
+ MakeSysError(TH1* tmp, TH1* cen, TH1* fwd, TF1* fit)
+ {
+ for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
+ Double_t tc = tmp->GetBinContent(i);
+ if (tc < 0.01) continue;
+ Double_t fc = fwd->GetBinContent(i);
+ Double_t cc = cen->GetBinContent(i);
+ Double_t sysErr = fFwdSysErr;
+ if (cc > .01 && fc > 0.01)
+ sysErr = (fFwdSysErr+fCenSysErr) / 2;
+ else if (cc > .01)
+ sysErr = fCenSysErr;
+ Double_t x = tmp->GetXaxis()->GetBinCenter(i);
+ Double_t y = fit->Eval(x);
+ tmp->SetBinContent(i, y);
+ tmp->SetBinError(i,sysErr*y);
+ }
+ TString name(tmp->GetName());
+ name.ReplaceAll("Merged", "SysError");
+ tmp->SetName(name);
+ tmp->SetFillColor(kGray);
+ tmp->SetFillStyle(3001);
+ tmp->SetMarkerStyle(0);
+ tmp->SetLineWidth(0);
+ }
+ void CorrectForward(TH1* h) const
+ {
+ if (!fRemoveOuters) return;
+
+ for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
+ Double_t eta = h->GetBinCenter(i);
+ if (TMath::Abs(eta) < 2.3) {
+ h->SetBinContent(i, 0);
+ h->SetBinError(i, 0);
+ }
+ }
+ }
+ void CorrectCentral(TH1* h) const
+ {
+ if (fClusterScale.IsNull()) return;
+ TString t(h->GetTitle());
+ Info("CorrectCentral", "Replacing Central with Tracklets in %s", t.Data());
+ t.ReplaceAll("Central", "Tracklets");
+ h->SetTitle(t);
+
+ TF1* cf = new TF1("clusterScale", fClusterScale);
+ for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
+ Double_t eta = h->GetBinCenter(i);
+ Double_t f = cf->Eval(eta);
+ Double_t c = h->GetBinContent(i);
+ if (f < .1) f = 1;
+ h->SetBinContent(i, c / f);
+ }
+ delete cf;
+ }
+
+
/* @} */
//__________________________________________________________________
- Bool_t fShowOthers; // Show other data
+ UShort_t fShowOthers; // Show other data
Bool_t fShowAlice; // Show ALICE published data
Bool_t fShowRatios; // Show ratios
Bool_t fShowLeftRight;// Show asymmetry
TH1* fTriggers; // Number of triggers
TH1* fTruth; // Pointer to truth
RangeParam* fRangeParam; // Parameter object for range zoom
-
+ Double_t fFwdSysErr; // Systematic error in forward range
+ Double_t fCenSysErr; // Systematic error in central range
+ Bool_t fRemoveOuters; // Whether to remove outers
+ TString fClusterScale; // Scaling of clusters to tracklets
};
+//____________________________________________________________________
+/**
+ * Function to calculate
+ * @f[
+ * g(x;A_1,A_2,\sigma_1,\sigma_2) =
+ * A_1\left(\frac{1}{2\pi\sigma_1}e^{(x/\sigma_1)^2} -
+ * A_2\frac{1}{2\pi\sigma_2}e^{(x/\sigma_2)^2}\right)
+ * @f]
+ *
+ * @param xp Pointer to x array
+ * @param pp Pointer to parameter array
+ *
+ * @return @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$
+ */
+Double_t myFunc(Double_t* xp, Double_t* pp)
+{
+ Double_t x = xp[0];
+ Double_t a1 = pp[0];
+ Double_t a2 = pp[1];
+ Double_t s1 = pp[2];
+ Double_t s2 = pp[3];
+ return a1*(TMath::Gaus(x, 0, s1) - a2 * TMath::Gaus(x, 0, s2));
+}
+
//=== Stuff for auto zooming =========================================
void UpdateRange(dNdetaDrawer::RangeParam* p)
{
*/
void
DrawdNdeta(const char* filename="forward_dndeta.root",
- Int_t flags=0xf,
const char* title="",
UShort_t rebin=5,
- Bool_t cutEdges=false,
+ UShort_t others=0x7,
+ UShort_t flags=0x7,
UShort_t sNN=0,
UShort_t sys=0,
UShort_t trg=0,
dNdetaDrawer* pd = new dNdetaDrawer;
dNdetaDrawer& d = *pd;
d.SetRebin(rebin);
- d.SetCutEdges(cutEdges);
d.SetTitle(title);
- d.SetShowOthers(flags & 0x1);
- d.SetShowAlice(flags & 0x2);
- d.SetShowRatios(flags & 0x4);
- d.SetShowLeftRight(flags & 0x8);
- d.SetShowRings(flags & 0x10);
+ d.SetShowOthers(others);
+ d.SetShowRatios(flags & 0x1);
+ d.SetShowLeftRight(flags & 0x2);
+ d.SetForwardSysError(flags & 0x4 ? 0.076 : 0);
+ d.SetShowRings(flags & 0x8);
+ d.SetCutEdges(flags & 0x10);
+ // d.fRemoveOuters = (flags & 0x20);
+ // d.fClusterScale = "1.06 -0.003*x +0.0119*x*x";
// Do the below if your input data does not contain these settings
if (sNN > 0) d.SetSNN(sNN); // Collision energy per nucleon pair (GeV)
if (sys > 0) d.SetSys(sys); // Collision system (1:pp, 2:PbPB)
*/
TGraphAsymmErrors* AliceCentralInelGt900()
{
+#if 0
// INEL>0 - p7741_d4x1y1 - Eur.Phys.J.C68:345-354,2010.
double x[] = { -0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3, 0.5, 0.7,
0.9 };
}
TGraphAsymmErrors* g = new TGraphAsymmErrors(np, x, y, exm, exp, eym, eyp);
+#else
+ // These are from JFGO
+ TGraphAsymmErrors *g = new TGraphAsymmErrors(15);
+ g->SetPoint(0,-1.5,4.12575);
+ g->SetPointError(0,0.1,0.1,0.0742967,0.0814571);
+ g->SetPoint(1,-1.3,3.91209);
+ g->SetPointError(1,0.1,0.1,0.0697701,0.0766199);
+ g->SetPoint(2,-1.1,3.98377);
+ g->SetPointError(2,0.1,0.1,0.0704503,0.0774795);
+ g->SetPoint(3,-0.9,4.00035);
+ g->SetPointError(3,0.1,0.1,0.0702388,0.0773433);
+ g->SetPoint(4,-0.7,3.87228);
+ g->SetPointError(4,0.1,0.1,0.067597,0.0745103);
+ g->SetPoint(5,-0.5,3.79613);
+ g->SetPointError(5,0.1,0.1,0.0659771,0.0727816);
+ g->SetPoint(6,-0.3,3.70489);
+ g->SetPointError(6,0.1,0.1,0.0642016,0.0708603);
+ g->SetPoint(7,-0.1,3.67423);
+ g->SetPointError(7,0.1,0.1,0.0635759,0.0701884);
+ g->SetPoint(8,0.1,3.72765);
+ g->SetPointError(8,0.1,0.1,0.0645004,0.071209);
+ g->SetPoint(9,0.3,3.72171);
+ g->SetPointError(9,0.1,0.1,0.064493,0.071182);
+ g->SetPoint(10,0.5,3.77428);
+ g->SetPointError(10,0.1,0.1,0.0655974,0.0723627);
+ g->SetPoint(11,0.7,3.91704);
+ g->SetPointError(11,0.1,0.1,0.0683783,0.0753716);
+ g->SetPoint(12,0.9,4.00674);
+ g->SetPointError(12,0.1,0.1,0.0703511,0.0774669);
+ g->SetPoint(13,1.1,3.97948);
+ g->SetPointError(13,0.1,0.1,0.0703744,0.077396);
+ g->SetPoint(14,1.3,3.99165);
+ g->SetPointError(14,0.1,0.1,0.0711888,0.078178);
+#endif
SetGraphAttributes(g, INELGt0, ALICE, false, "alice_inelgt900",
"ALICE INEL>0 (publ.)");
return g;
-
}
//____________________________________________________________________
*/
TGraphAsymmErrors* AliceCentralInelGt2360()
{
+#if 0
// INEL>0 - p7741_d5x1y1 - Eur.Phys.J.C68:345-354,2010.
double x[] = { -0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3, 0.5, 0.7, 0.9 };
double exm[] = { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1 };
}
TGraphAsymmErrors* g = new TGraphAsymmErrors(np, x, y, exm, exp, eym, eyp);
+#else
+ // These are from JFGO
+ TGraphAsymmErrors* g = new TGraphAsymmErrors(15);
+ g->SetPoint(0,-1.5,4.79047);
+ g->SetPointError(0,0.1,0.1,0.0844278,0.109947);
+ g->SetPoint(1,-1.3,4.91068);
+ g->SetPointError(1,0.1,0.1,0.0856751,0.112038);
+ g->SetPoint(2,-1.1,4.87386);
+ g->SetPointError(2,0.1,0.1,0.0842846,0.110628);
+ g->SetPoint(3,-0.9,4.91365);
+ g->SetPointError(3,0.1,0.1,0.084339,0.111049);
+ g->SetPoint(4,-0.7,4.7601);
+ g->SetPointError(4,0.1,0.1,0.0812087,0.107203);
+ g->SetPoint(5,-0.5,4.63355);
+ g->SetPointError(5,0.1,0.1,0.078687,0.104079);
+ g->SetPoint(6,-0.3,4.63885);
+ g->SetPointError(6,0.1,0.1,0.0785337,0.104014);
+ g->SetPoint(7,-0.1,4.55439);
+ g->SetPointError(7,0.1,0.1,0.0769842,0.10203);
+ g->SetPoint(8,0.1,4.55087);
+ g->SetPointError(8,0.1,0.1,0.0769246,0.101951);
+ g->SetPoint(9,0.3,4.64118);
+ g->SetPointError(9,0.1,0.1,0.0785732,0.104066);
+ g->SetPoint(10,0.5,4.66172);
+ g->SetPointError(10,0.1,0.1,0.0791652,0.104711);
+ g->SetPoint(11,0.7,4.81871);
+ g->SetPointError(11,0.1,0.1,0.0822086,0.108523);
+ g->SetPoint(12,0.9,4.88193);
+ g->SetPointError(12,0.1,0.1,0.0837944,0.110332);
+ g->SetPoint(13,1.1,4.89068);
+ g->SetPointError(13,0.1,0.1,0.0845754,0.111009);
+ g->SetPoint(14,1.3,5.05663);
+ g->SetPointError(14,0.1,0.1,0.0882216,0.115368);
+#endif
+
SetGraphAttributes(g, INELGt0, ALICE, false, "alice_inelgt2360",
"ALICE INEL>0 (publ.)");
return g;
*/
TGraphAsymmErrors* AliceCentralInelGt7000()
{
+#if 0
// INEL>0 - p7741_d6x1y1 - Eur.Phys.J.C68:345-354,2010.
// Plot: p7741_d6x1y1
double x[] = { -0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3, 0.5, 0.7, 0.9 };
}
TGraphAsymmErrors* g = new TGraphAsymmErrors(np, x, y, exm, exp, eym, eyp);
+#else
+ // These are from JFGO
+ TGraphAsymmErrors* g = new TGraphAsymmErrors(15);
+ g->SetPoint(0,-1.5,6.28573);
+ g->SetPointError(0,0.1,0.1,0.125928,0.215392);
+ g->SetPoint(1,-1.3,6.25573);
+ g->SetPointError(1,0.1,0.1,0.124352,0.213795);
+ g->SetPoint(2,-1.1,6.28779);
+ g->SetPointError(2,0.1,0.1,0.124143,0.214399);
+ g->SetPoint(3,-0.9,6.21881);
+ g->SetPointError(3,0.1,0.1,0.122079,0.211642);
+ g->SetPoint(4,-0.7,6.0728);
+ g->SetPointError(4,0.1,0.1,0.118661,0.206355);
+ g->SetPoint(5,-0.5,6.011);
+ g->SetPointError(5,0.1,0.1,0.117043,0.204019);
+ g->SetPoint(6,-0.3,5.84071);
+ g->SetPointError(6,0.1,0.1,0.11346,0.198086);
+ g->SetPoint(7,-0.1,5.8532);
+ g->SetPointError(7,0.1,0.1,0.113569,0.198433);
+ g->SetPoint(8,0.1,5.84811);
+ g->SetPointError(8,0.1,0.1,0.11347,0.198261);
+ g->SetPoint(9,0.3,5.91022);
+ g->SetPointError(9,0.1,0.1,0.11481,0.200444);
+ g->SetPoint(10,0.5,6.00649);
+ g->SetPointError(10,0.1,0.1,0.116955,0.203866);
+ g->SetPoint(11,0.7,6.17115);
+ g->SetPointError(11,0.1,0.1,0.120583,0.209697);
+ g->SetPoint(12,0.9,6.2645);
+ g->SetPointError(12,0.1,0.1,0.122976,0.213197);
+ g->SetPoint(13,1.1,6.36448);
+ g->SetPointError(13,0.1,0.1,0.125657,0.217014);
+ g->SetPoint(14,1.3,6.39489);
+ g->SetPointError(14,0.1,0.1,0.127118,0.218551);
+#endif
SetGraphAttributes(g, INELGt0, ALICE, false, "alice_inelgt7000",
"ALICE INEL>0 (publ.)");
return g;
return g;
}
+//____________________________________________________________________
+TGraphAsymmErrors*
+GetSingle(UShort_t which,
+ UShort_t sys,
+ UShort_t energy,
+ UShort_t type=0x1,
+ UShort_t centLow=0,
+ UShort_t centHigh=0)
+{
+ TGraphAsymmErrors* ret = 0;
+ if (sys == 1) {
+ if (TMath::Abs(energy-900) < 10) {
+ switch (type) {
+ case 1: // INEL
+ switch (which) {
+ case PYTHIA: ret = Pythia900INEL(); break;
+ case UA5: ret = UA5Inel(false); break;
+ case UA5+10: ret = UA5Inel(true); break;
+ case ALICE: ret = AliceCentralInel900(); break;
+ }
+ break;
+ case 2: // INEL>0
+ switch (which) {
+ case ALICE: ret = AliceCentralInelGt900(); break;
+ }
+ break;
+ case 4: // NSD
+ switch (which) {
+ case PYTHIA: ret = Pythia900NSD(); break;
+ case UA5: ret = UA5Nsd(false); break;
+ case UA5+10: ret = UA5Nsd(true); break;
+ case ALICE: ret = AliceCentralNsd900(); break;
+ case CMS: ret = CMSNsd900(); break;
+ }
+ break;
+ } // type
+ }
+ else if (TMath::Abs(energy-2360) < 10) {
+ switch (type) {
+ case 1: // INEL
+ switch (which) {
+ case ALICE: ret = AliceCentralInel2360(); break;
+ }
+ break;
+ case 2: // INEL > 0
+ switch (which) {
+ case ALICE: ret = AliceCentralInelGt2360(); break;
+ }
+ break;
+ case 4: // NSD
+ switch (which) {
+ case ALICE: ret = AliceCentralNsd2360(); break;
+ case CMS: ret = CMSNsd2360(); break;
+ }
+ break;
+ }
+ }
+ else if (TMath::Abs(energy-7000) < 10) {
+ switch (type) {
+ case 1: ret = 0; break;
+ case 2: // INEL > 0
+ switch (which) {
+ case ALICE: ret = AliceCentralInelGt7000(); break;
+ }
+ break;
+ case 4: // NSD
+ switch (which) {
+ case CMS: ret = CMSNsd7000(); break;
+ }
+ break;
+ }
+ }
+ }
+ return ret;
+}
+
+
//____________________________________________________________________
/**
* Get a multi graph of data for a given energy and trigger type
UShort_t type=0x1,
UShort_t centLow=0,
UShort_t centHigh=0,
- bool aliceOnly=false,
- bool nomc=false)
+ UShort_t which=0x7)
{
TMultiGraph* mp = new TMultiGraph(Form("dndeta_%dGeV_%d_%03d_%03d",
energy, type, centLow, centHigh),"");
TString en;
TString sn;
TString cn;
- if (sys == 1) {
- sn = ", pp(p#bar{p})";
- if (energy < 1000)
- en = Form(", #sqrt{s}=%dGeV", energy);
+ bool ua5 = (which & (1 << UA5)); // 0x1
+ bool cms = (which & (1 << CMS)); // 0x2
+ bool alice = (which & (1 << ALICE)); // 0x4
+ bool pythia = (which & (1 << PYTHIA)); // 0x8
+
+ en.Append(Form(", #sqrt{s%s}=", sys == 1 ? "" : "_{NN}"));
+ if (energy < 1000)
+ en.Append(Form("%dGeV", energy));
+ else {
+ if (energy % 1000 == 0)
+ en.Append(Form("%dTeV", energy/1000));
else
- en = Form(", #sqrt{s}=%f4.2TeV", float(energy)/1000);
+ en.Append(Form("%4.2fTeV", float(energy)/1000));
+ }
+
+ if (sys == 1) {
if (!(type & 0x7))
Warning("GetData", "Unknown trigger mask 0x%x", type);
- if (TMath::Abs(energy-900) < 10) {
- if (type & 0x1) {
- tn.Append(" INEL");
- if (!aliceOnly && !nomc) mp->Add(Pythia900INEL());
- if (!aliceOnly) mp->Add(UA5Inel(false));
- if (!aliceOnly) mp->Add(UA5Inel(true));
- mp->Add(AliceCentralInel900());
- }
- if (type & 0x4) {
- tn.Append(" NSD");
- if (!aliceOnly && !nomc) mp->Add(Pythia900NSD());
- if (!aliceOnly) mp->Add(UA5Nsd(false));
- if (!aliceOnly) mp->Add(UA5Nsd(true));
- mp->Add(AliceCentralNsd900());
- if (!aliceOnly) mp->Add(CMSNsd900());
- }
- if (type & 0x2) {
- tn.Append(" INEL>0");
- mp->Add(AliceCentralInelGt900());
- }
- }
- else if (TMath::Abs(energy-2360) < 10) {
- if (type & 0x1) {
- tn.Append(" INEL");
- mp->Add(AliceCentralInel2360());
- }
- if (type & 0x4) {
- tn.Append(" NSD");
- mp->Add(AliceCentralNsd2360());
- if (!aliceOnly) mp->Add(CMSNsd2360());
- }
- if (type & 0x2) {
- tn.Append(" INEL>0");
- mp->Add(AliceCentralInelGt2360());
- }
- }
- else if (TMath::Abs(energy-7000) < 10) {
- if (type & 0x1) {
- tn.Append(" INEL");
- }
- if (type & 0x4) {
- tn.Append(" NSD");
- if (!aliceOnly) mp->Add(CMSNsd7000());
- }
- if (type & 0x2) {
- tn.Append(" INEL>0");
- mp->Add(AliceCentralInelGt7000());
- }
- }
-#if 0
- else
+ if (!(TMath::Abs(energy-900) < 10 ||
+ TMath::Abs(energy-2360) < 10 ||
+ TMath::Abs(energy-7000) < 10)) {
Warning("GetData", "No other results for sys=%d, energy=%d",
sys, energy);
-#endif
+ return 0;
+ }
+
+ sn = "pp";
+
+ if (type & 0x1) tn.Append("INEL");
+ if (type & 0x2) { if (!tn.IsNull()) tn.Append("|"); tn.Append("INEL>0"); }
+ if (type & 0x4) { if (!tn.IsNull()) tn.Append("|"); tn.Append("NSD"); }
+
+ Bool_t seenUA5 = false;
+ for (Int_t i = 0; i < 3; i++) {
+ UShort_t mask = (1 << i);
+ if ((type & mask) == 0) continue;
+ TGraphAsymmErrors* gUAp = (ua5 ? GetSingle(UA5, sys,energy,mask): 0);
+ TGraphAsymmErrors* gUAn = (ua5 ? GetSingle(UA5+10,sys,energy,mask): 0);
+ TGraphAsymmErrors* gCMS = (cms ? GetSingle(CMS, sys,energy,mask): 0);
+ TGraphAsymmErrors* gALI = (alice ? GetSingle(ALICE, sys,energy,mask): 0);
+ TGraphAsymmErrors* gPYT = (pythia? GetSingle(PYTHIA,sys,energy,mask): 0);
+ if (gUAp) mp->Add(gUAp);
+ if (gUAn) mp->Add(gUAn);
+ if (gCMS) mp->Add(gCMS);
+ if (gALI) mp->Add(gALI);
+ if (gPYT) mp->Add(gPYT);
+ if (gUAp || gUAn) seenUA5 = true;
+ }
+ if (seenUA5) sn.Append("(p#bar{p})");
}
else if (sys == 2) {
// Nothing for PbPb so far
cn = Form(", %d%%-%d%% central", centLow, centHigh);
- sn = ", PbPb";
- if (energy < 1000)
- en = Form(", #sqrt{s_{NN}}=%dGeV", energy);
- else
- en = Form(", #sqrt{s_{NN}}=%f4.2TeV", float(energy)/1000);
+ sn = "PbPb";
// Warning("GetData", "No other data for PbPb yet");
}
else
Warning("GetData", "Unknown system %d", sys);
- TString tit(Form("1/N dN_{ch}/d#eta%s%s%s%s",
- sn.Data(), en.Data(), tn.Data(), cn.Data()));
- mp->SetTitle(tit.Data());
+
if (!mp->GetListOfGraphs() || mp->GetListOfGraphs()->GetEntries() <= 0) {
delete mp;
mp = 0;
+ return 0;
}
+ TString tit(Form("%s%s, %s%s",
+ sn.Data(), en.Data(), tn.Data(), cn.Data()));
+ mp->SetTitle(tit.Data());
return mp;
}
* - 0x4 NSD
* @param centLow Low centrality cut (only for PbPB)
* @param centHigh High centrality cut (only for PbPB)
- * @param aliceOnly Only return other ALICE data
+ * @param alice Only return other ALICE data
*
* @ingroup pwg2_forward_otherdata
*/
void
OtherData(UShort_t sys=1,
- UShort_t energy=900,
- UShort_t type=0x1,
- UShort_t centLow=0,
- UShort_t centHigh=5,
- bool aliceOnly=false)
+ UShort_t energy=900,
+ UShort_t type=0x1,
+ UShort_t centLow=0,
+ UShort_t centHigh=5,
+ UShort_t which=0x7)
{
- TMultiGraph* mp = GetData(sys, energy, type, centLow, centHigh, aliceOnly);
+ TMultiGraph* mp = GetData(sys, energy, type, centLow, centHigh, which);
if (!mp) return;
gStyle->SetTitleX(0.1);
gStyle->SetTitleFillColor(kBlack);
gStyle->SetTitleFontSize(0.02);
- gStyle->SetOptTitle(1);
+ gStyle->SetOptTitle(0);
gStyle->SetOptStat(0);
TCanvas* c = new TCanvas("c", "dN/deta", 800, 600);
c->SetFillColor(0);
c->SetBorderSize(0);
c->SetBorderMode(0);
- c->SetRightMargin(0.05);
- c->SetTopMargin(0.05);
+ c->SetRightMargin(0.02);
+ c->SetTopMargin(0.02);
mp->SetMinimum(0);
mp->Draw("ap");
- if (mp->GetXaxis())
+ if (mp->GetXaxis()) {
mp->GetXaxis()->SetTitle("#eta");
- if (mp->GetYaxis())
- mp->GetYaxis()->SetTitle("#frac{1}{N} #frac{dN_{ch}}{#eta}");
-
- TLegend* l = c->BuildLegend(0.3, 0.15, 0.7, 0.5);
+ mp->GetXaxis()->SetTitleFont(12);
+ mp->GetXaxis()->SetLabelFont(132);
+ }
+ if (mp->GetYaxis()) {
+ mp->GetYaxis()->SetTitle("#frac{1}{N} #frac{dN_{#font[132]{ch}}}{d#eta}");
+ mp->GetYaxis()->SetTitleFont(12);
+ mp->GetYaxis()->SetLabelFont(132);
+ }
+ TLegend* l = c->BuildLegend(0.3, 0.15, 0.7, 0.5,
+ mp->GetTitle());
l->SetFillColor(0);
+ l->SetTextFont(132);
l->SetBorderSize(0);
-
+ TLegendEntry* h = static_cast<TLegendEntry*>(l->GetListOfPrimitives()->At(0));
+ if (h) {
+ h->SetTextFont(22);
+ }
c->cd();
}
#!/bin/bash
+# General options
+batch=0
+proof=0
+cent=0
+dopass1=0
+dopass2=0
+dopass3=0
+max_rotate=10
+prog=aliroot
+
+# Input (data and scripts)
ana=$ALICE_ROOT/PWG2/FORWARD/analysis2
esddir="."
nev=-1
-rebin=1
+
+# Pass 1 options
+mc=0
+scheme="full"
+
+# Pass 2 event selection
vzmin=-10
vzmax=10
-batch=0
-proof=0
-mc=0
type=INEL
-cms=900
-hhd=1
-comp=1
-cent=0
+mcfilename="none"
+
+# Pass 3 visualisation options
+rebin=1
tit=
-others=0
-published=1
+others=
ratios=1
asymm=1
-scheme="full"
-dopass1=0
-dopass2=0
-dopass3=0
+rings=0
+syserr=1
+
+# Derived variables
pass1=MakeAOD.C
pass2=MakedNdeta.C
-pass3=DrawdNdeta.C++g
+pass3=DrawdNdeta.C++
output1=forward.root
output2=forward_dndeta.root
outputs1="${output1} AliAOD.root event_stat.root EventStat_temp.root"
outputs2="${output2}"
gdb_script=$ALICE_ROOT/PWG2/FORWARD/analysis2/gdb_cmds
-max_rotate=10
name=`date +analysis%Y%m%d_%H%M`
pass2dir=./
-mcfilename="none"
#_____________________________________________________________________
# Print usage
Do Pass1 and Pass2 on ESD files in current directory.
Options:
+ General options:
-h,--help This help
- -n,--events N Number of events ($nev)
- -1,--pass1 Run pass 1, only AOD ($dopass1)
- -2,--pass2 Run pass 2, only Hists ($dopass2)
- -3,--pass3 Draw results ($dopass3)
- -v,--vz-min CM Minimum value of vz ($vzmin)
- -V,--vz-max CM Maximum value of vz ($vzmax)
- -t,--trigger TYPE Select trigger TYPE ($type)
-b,--batch Do batch processing ($batch)
-P,--proof NWORKERS Run in PROOF(Lite) mode ($proof)
- -M,--mc Run over MC data ($mc)
+ -1,--pass1 Run pass 1 (AOD) ($dopass1)
+ -2,--pass2 Run pass 2 (Hists) ($dopass2)
+ -3,--pass3 Run pass 3 (Vizualisation) ($dopass3)
+ -n,--events N Number of events ($nev)
-g,--gdb Run in GDB mode ($gdb)
+ -N,--name STRING Name of analysis ($name)
-E,--eloss Run energy loss script
- -r,--rebin Rebin factor ($rebin)
+ -A,--script-dir Script directory ($ana)
+ -a,--program Program to use ($prog)
+ Pass 1 options and event selection
-C,--use-centrality Run centrality task ($cent)
- -O,--show-older Show older data ($others)
- -J,--show-published Show ALICE published data ($published)
+ -M,--mc Run over MC data ($mc)
+ -I,--input-dir PATH Path to input ESD data ($esddir)
+ Pass 2 options
+ -v,--vz-min CM Minimum value of vz ($vzmin)
+ -V,--vz-max CM Maximum value of vz ($vzmax)
+ -t,--trigger TYPE Select trigger TYPE ($type)
+ -S,--scheme SCHEME Normalisation scheme ($scheme)
+ -q,--mcfilename Result of MC analysis ($mcfilename)
+ Pass 3 options
+ -r,--rebin Rebin factor ($rebin)
+ -O,--others WHICH Show other data ($others)
+ -J,--show-rings Show individual rings ($rings)
-R,--show-ratios Show ratios to other data ($ratios)
-Z,--show-asymmetry Show asymmetry ($asymm)
- -S,--scheme SCHEME Normalisation scheme ($scheme)
+ -G,--show-syserror Show systematic errors ($syserr)
-T,--title STRING Title on plots ($tit)
- -q,--mcfilename MC dndeta filename ($mcfilename)
- -N,--name STRING Name of analysis ($name)
- -I,--input-dir PATH Path to input ESD data ($esddir)
TYPE is a comma or space separated list of
TRIGGER Trigger efficiency
FULL Same as EVENTLEVEL,BACKGROUND,SHAPE,TRIGGER
+WHICH is a comma or space separated list of
+
+ UA5 UA5 data (900GeV, ppbar, INEL, NSD)
+ CMS CMS data
+ ALICE ALICE published data
+ PYTHIA Event generator data
+
If NWORKERS is 0, then the analysis will be run in local mode.
EOF
}
EOF
# --- Run AliROOT ------------------------------------------------
- echo "Running aliroot $opts ${script}${args}"
+ echo "Running ${prog} $opts ${script}${args}"
if test $isbatch -gt 0 || test $notLast -gt 0 ; then
- aliroot ${opts} ${script}"${args}" 2>&1 | tee ${log}
+ ${prog} ${opts} ${script}"${args}" 2>&1 | tee ${log}
else
- aliroot ${opts} ${script}"${args}"
+ ${prog} ${opts} ${script}"${args}"
fi
fail=$?
# Loop over arguments
while test $# -gt 0 ; do
case $1 in
- -h|--help) usage ; exit 0;;
- -n|--events) nev=$2 ; shift ;;
- -3|--pass3|-D|--draw) dopass3=`toggle $dopass3` ;;
- -2|--pass2|-H|--hist) dopass2=`toggle $dopass2` ;;
- -1|--pass1|-A|--aod) dopass1=`toggle $dopass1` ;;
+ # General options
+ -h|--help) usage ; exit 0 ;;
-b|--batch) batch=`toggle $batch` ;;
- -P|--proof) proof=$2 ; shift ;;
- -C|--use-centrality) cent=`toggle $cent` ;;
- -M|--mc) mc=`toggle $mc` ;;
- -g|--gdb) gdb=`toggle $gdb` ;;
- -r|--rebin) rebin=$2 ; shift ;;
- -v|--vz-min) vzmin=$2 ; shift ;;
- -V|--vz-max) vzmax=$2 ; shift ;;
+ -P|--proof) proof=$2 ; shift ;;
+ -1|--pass1|--aod) dopass1=`toggle $dopass1` ;;
+ -2|--pass2|--hist) dopass2=`toggle $dopass2` ;;
+ -3|--pass3|--draw) dopass3=`toggle $dopass3` ;;
+ -n|--events) nev=$2 ; shift ;;
+ -g|--gdb) gdb=`toggle $gdb` ;;
+ -N|--name) name=$2 ; shift ;;
+ -A|--script-dir) ana=$2; shift ;;
+ -a|--program) prog=$2 ; shift ;;
-E|--eloss) pass1=MakeELossFits.C
pass2=scripts/ExtractELoss.C
pass3=scripts/DrawAnaELoss.C
outputs2=""
dopass2=1
;;
- -O|--show-older) others=`toggle $others` ;;
- -J|--show-published) published=`toggle $published` ;;
+ # Pass 1 options
+ -C|--use-centrality) cent=`toggle $cent` ;;
+ -M|--mc) mc=`toggle $mc` ;;
+ -I|--input-dir) esddir=$2 ; shift ;;
+ # Pass 2 options
+ -v|--vz-min) vzmin=$2 ; shift ;;
+ -V|--vz-max) vzmax=$2 ; shift ;;
+ -S|--scheme) scheme=`echo $2 | tr ' ' ','` ; shift ;;
+ -t|--type) type=$2 ; shift ;;
+ -q|--mcname) mcfilename=$2 ; shift ;;
+ # PAss 3 options
+ -r|--rebin) rebin=$2 ; shift ;;
+ -O|--others) others=`echo $2 | tr ',' ' '` ; shift ;;
-R|--show-ratios) ratios=`toggle $ratios` ;;
+ -J|--show-rings) rings=`toggle $rings` ;;
-Z|--show-asymmetry) asymm=`toggle $asymm` ;;
- -S|--scheme) scheme=`echo $2 | tr ' ' ','` ; shift ;;
+ -G|--show-syserror) syserr=`toggle $syserr` ;;
-T|--title) tit=$2 ; shift ;;
- -q|--mcname) mcfilename=$2 ; shift ;;
- -N|--name) name=$2 ; shift ;;
- -I|--input-dir) esddir=$2 ; shift ;;
- -t|--type)
- #if test "x$type" = "x" ; then type=$2 ; else type="$type|$2"; fi
- type=$2
- shift ;;
*) echo "$0: Unknown option '$1'" >> /dev/stderr ; exit 1 ;;
esac
shift
if test $dopass3 -gt 0 ; then
tit=`echo $tit | tr ' ' '@'`
flags=0
- if test $others -gt 0 ; then let flags=$(($flags|0x1)); fi
- if test $published -gt 0 ; then let flags=$(($flags|0x2)); fi
- if test $ratios -gt 0 ; then let flags=$(($flags|0x4)); fi
- if test $asymm -gt 0 ; then let flags=$(($flags|0x8)); fi
-
- echo "Running with Title=$tit"
+ if test $ratios -gt 0 ; then let flags=$(($flags|0x1)); fi
+ if test $asymm -gt 0 ; then let flags=$(($flags|0x2)); fi
+ if test $syserr -gt 0 ; then let flags=$(($flags|0x4)); fi
+ if test $rings -gt 0 ; then let flags=$(($flags|0x8)); fi
+ which=0
+ others=`echo $others | tr ',' ' ' | tr '[a-z]' '[A-Z]'`
+ for i in $others ; do
+ case $i in
+ UA5) let which=$(($which|0x1)) ;;
+ CMS) let which=$(($which|0x2)) ;;
+ ALICE) let which=$(($which|0x4)) ;;
+ PYTHIA|MC) let which=$(($which|0x8)) ;;
+ esac
+ done
- args="(\"${pass2dir}${output2}\",${flags},\"$tit\",$rebin)"
+ args="(\"${pass2dir}${output2}\",\"$tit\",$rebin,${which},${flags})"
if test "x$pass1" = "xMakeELossFits.C" ; then
args="(\"${pass2dir}${output1}\")"
fi
+
Double_t
-DrawRingELossPoisson(TList* p, UShort_t d, Char_t r)
+DrawRingELossPoisson(TList* p, UShort_t d, Char_t r,
+ Double_t xmin, Double_t xmax)
{
if (!p) return;
gPad->SetGridx();
gPad->SetLogz();
gPad->SetFillColor(0);
+ if (xmax < 0) xmax = corr->GetXaxis()->GetXmax();
+ corr->GetXaxis()->SetRangeUser(xmin,xmax);
+ corr->GetYaxis()->SetRangeUser(xmin,xmax);
corr->SetTitle(Form("FMD%d%c",d,r));
corr->Draw("colz");
- TLine* l = new TLine(0,0,corr->GetXaxis()->GetXmax(),
- corr->GetYaxis()->GetXmax());
+ // Calculate the linear regression
+ Double_t dx = (xmax-xmin);
+ Double_t rxy = corr->GetCorrelationFactor();
+ Double_t sx = corr->GetRMS(1);
+ Double_t sy = corr->GetRMS(2);
+ Double_t sx2 = sx*sx;
+ Double_t sy2 = sy*sy;
+ Double_t cxy = corr->GetCovariance();
+ Double_t mx = corr->GetMean(1);
+ Double_t my = corr->GetMean(2);
+ Double_t delta = 1;
+#if 0
+ Double_t beta = rxy * sy / sx;
+#else
+ Double_t beta = ((sy2 - delta*sx2 +
+ TMath::Sqrt(TMath::Power(sy2-delta*sx2,2) +
+ 4*delta*cxy*cxy))
+ / 2 / cxy);
+#endif
+ Double_t alpha = my - beta * mx;
+ TF1* f = new TF1("f", "[0]+[1]*x", xmin, xmax);
+ f->SetParameters(alpha, beta);
+ f->SetLineWidth(1);
+ f->SetLineStyle(1);
+ f->Draw("same");
+
+ TLine* l = new TLine(xmin,xmin,xmax,xmax);
l->SetLineWidth(1);
+ l->SetLineStyle(2);
l->SetLineColor(kBlack);
l->Draw();
// corr->GetXaxis()->SetRangeUser(0, 2);
- Info("", "FMD%d%c correlation coefficient: %9.5f%%",
- d, r, 100*corr->GetCorrelationFactor());
+ Info("", "FMD%d%c correlation coefficient: %9.5f%% "
+ "line y = %f + %f * x", d, r, 100*rxy, alpha, beta);
+
+ Double_t x = xmin + dx * .1;
+ Double_t y = xmax - dx * .2;
+ TLatex* ltx = new TLatex(x, y, "Deming regression: y=#alpha+#beta x");
+ ltx->SetTextAlign(13);
+ ltx->SetTextSize(0.06);
+ // ltx->SetNDC();
+ ltx->Draw();
+ ltx->SetTextSize(0.05);
+ x = xmin + dx / 8;
+ y = xmax - dx * .25;
+ ltx->DrawLatex(x, y, Form("#alpha=%5.3f", alpha));
+ y = xmax - dx * .3;
+ ltx->DrawLatex(x, y, Form("#beta= %5.3f", beta));
+
+ TLinearFitter* fitter = new TLinearFitter(1);
+ fitter->SetFormula("1 ++ x");
+ for (Int_t i = 1; i <= corr->GetNbinsX(); i++) {
+ x = corr->GetXaxis()->GetBinCenter(i);
+ if (x < xmin || x > xmax) continue;
+ for (Int_t j = 1; j <= corr->GetNbinsY(); j++) {
+ y = corr->GetYaxis()->GetBinCenter(j);
+ if (y < xmin || y > xmax) continue;
+ Double_t c = corr->GetBinContent(i,j);
+ if (c < .1) continue;
+ fitter->AddPoint(&x, y, c);
+ }
+ }
+ Bool_t robust = false;
+ if (robust)
+ fitter->EvalRobust();
+ else
+ fitter->Eval();
+ for (Int_t i = 0; i < 2; i++) {
+ std::cout << i << "\t"
+ << fitter->GetParName(i) << "\t"
+ << fitter->GetParameter(i) << "\t";
+ if (!robust)
+ std::cout << fitter->GetParError(i);
+ std::cout << std::endl;
+ }
+ Double_t chi2 = fitter->GetChisquare();
+ Int_t ndf = (fitter->GetNpoints() -
+ fitter->GetNumberFreeParameters() );
+ std::cout << chi2 << '/' << ndf << '=' << chi2 / ndf << std::endl;
+
gPad->cd();
return corr->GetCorrelationFactor();
}
void
-DrawELossPoisson(const char* filename="forward.root")
+DrawELossPoisson(const char* filename="forward.root",
+ Double_t xmax=-1,
+ Double_t xmin=0)
{
gStyle->SetPalette(1);
gStyle->SetOptFit(0);
c->Divide(3, 2, 0, 0);
Double_t corrs[5];
- c->cd(1); corrs[0] = DrawRingELossPoisson(dc, 1, 'I');
- c->cd(2); corrs[1] = DrawRingELossPoisson(dc, 2, 'I');
- c->cd(5); corrs[2] = DrawRingELossPoisson(dc, 2, 'O');
- c->cd(3); corrs[3] = DrawRingELossPoisson(dc, 3, 'I');
- c->cd(6); corrs[4] = DrawRingELossPoisson(dc, 3, 'O');
+ c->cd(1); corrs[0] = DrawRingELossPoisson(dc, 1, 'I', xmin, xmax);
+ c->cd(2); corrs[1] = DrawRingELossPoisson(dc, 2, 'I', xmin, xmax);
+ c->cd(5); corrs[2] = DrawRingELossPoisson(dc, 2, 'O', xmin, xmax);
+ c->cd(3); corrs[3] = DrawRingELossPoisson(dc, 3, 'I', xmin, xmax);
+ c->cd(6); corrs[4] = DrawRingELossPoisson(dc, 3, 'O', xmin, xmax);
TVirtualPad* p = c->cd(4);
p->SetTopMargin(0.05);
hc->SetFillColor(kRed+1);
hc->SetFillStyle(3001);
hc->SetMinimum(0.0);
- hc->SetMaximum(1.1);
+ hc->SetMaximum(1.3);
hc->GetXaxis()->SetBinLabel(1,"FMD1i"); hc->SetBinContent(1,corrs[0]);
hc->GetXaxis()->SetBinLabel(2,"FMD2i"); hc->SetBinContent(2,corrs[1]);
hc->GetXaxis()->SetBinLabel(3,"FMD2o"); hc->SetBinContent(3,corrs[2]);
// TH2D* highCuts = static_cast<TH2D*>(dc->FindObject("highCuts"));
// if (highCuts) highCuts->Draw("colz");
c->cd();
-
+ c->SaveAs("elossVsPoisson.png");
}
--- /dev/null
+void
+DrawRubensCorr(const char* fname="rubensRatio.root",
+ const char* hname = "dNdEtaCor1D_cls")
+{
+ TFile* f = TFile::Open(fname, "READ");
+ if (!f) {
+ Error("DrawRubensCorr", "%s not found", fname);
+ return;
+ }
+ TH1D* h = static_cast<TH1D*>(f->Get(hname));
+ if (!h) {
+ Error("DrawRubensCorr", "%s not found in %s", hname, fname);
+ return;
+ }
+
+ gStyle->SetOptTitle(0);
+ gStyle->SetOptStat(0);
+
+ TCanvas* c = new TCanvas("c", "c", 800, 800);
+ c->SetFillColor(0);
+ c->SetFillStyle(0);
+ c->SetBorderSize(0);
+ c->SetBorderMode(0);
+ c->SetRightMargin(0.03);
+ c->SetLeftMargin(0.15);
+ c->SetTopMargin(0.03);
+
+ h->GetXaxis()->SetTitleFont(132);
+ h->GetXaxis()->SetLabelFont(132);
+ h->GetXaxis()->SetDecimals();
+ h->GetYaxis()->SetTitleFont(132);
+ h->GetYaxis()->SetLabelFont(132);
+ h->GetYaxis()->SetDecimals();
+ h->GetYaxis()->SetTitleOffset(1.8);
+ h->SetYTitle("Clusters/Tracklets");
+ h->Draw();
+
+ TF1* g = new TF1("g", "pol2", -2.2, 2.2);
+ g->SetLineWidth(2);
+ h->Fit(g, "NQ");
+ g->Draw("same");
+ Info("DrawRubensCorr", "Result of Fit:\n"
+ " chi^2/NDF = %7.4f / %2d = %7.4f\n"
+ " A = %+4.2f +/- %5.3f\n"
+ " B = %+5.3f +/- %6.4f\n"
+ " C = %+6.4f +/- %7.5f\n"
+ " f(x)=%4.2f%+5.3f*x%+6.4f*x*x",
+ g->GetChisquare(), g->GetNDF(), g->GetChisquare()/g->GetNDF(),
+ g->GetParameter(0), g->GetParError(0),
+ g->GetParameter(1), g->GetParError(1),
+ g->GetParameter(2), g->GetParError(2),
+ g->GetParameter(0), g->GetParameter(1), g->GetParameter(2));
+
+ TLatex* ltx = new TLatex(.2, .9, Form("f(x)=%4.2f%+5.3fx%+6.4fx^{2}",
+ g->GetParameter(0),
+ g->GetParameter(1),
+ g->GetParameter(2)));
+ ltx->SetNDC();
+ ltx->SetTextFont(132);
+ ltx->Draw();
+
+ c->SaveAs("rubens_corr.png");
+}
+
--- /dev/null
+THStack*
+GetStack(const TList& forward, const char* sub, const char* name)
+{
+ TList* lsub = &forward;
+
+ if (sub && sub[0] != '\0')
+ lsub = static_cast<TList*>(forward.FindObject(sub));
+
+ if (!lsub) {
+ Warning("GetStack", "Sub list %s not found in %s", sub, forward.GetName());
+ return 0;
+ }
+ THStack* ret = static_cast<THStack*>(lsub->FindObject(name));
+ if (!ret)
+ Warning("GetStack" "Stack %s not found in %s", name, sub);
+ return ret;
+}
+
+TH1*
+Rebin(TH1* h, Int_t rebin)
+{
+ if (rebin <= 1) return h;
+ h->Rebin(rebin);
+ h->Scale(1. / rebin);
+ return h;
+}
+
+TH1*
+Ratio(const TH1* h1, const TH1* h2)
+{
+ if (!h1) return;
+ if (!h2) return;
+
+ TH1* copy = static_cast<TH1*>(h2->Clone("tmp"));
+ copy->SetName(Form("%s_%s", h2->GetName(), h1->GetName()));
+ copy->SetTitle(Form("%s/%s", h2->GetTitle(), h1->GetTitle()));
+ copy->SetDirectory(0);
+ copy->Divide(h1);
+
+ return copy;
+}
+
+Int_t
+Ratio(THStack* r, const THStack* h1, const THStack* h2)
+{
+ if (!h1) return 0;
+ if (!h2) return 0;
+
+ int n1 = h1->GetHists()->GetEntries();
+ int n2 = h2->GetHists()->GetEntries();
+ int nH = 0;
+ for (int i = 0; i < n1 && i < n2; i++) {
+ TH1* hh1 = static_cast<TH1*>(h1->GetHists()->At(i));
+ TH1* hh2 = static_cast<TH1*>(h2->GetHists()->At(i));
+ TH1* h = Ratio(hh1, hh2);
+ if (!h) continue;
+ nH++;
+ r->Add(h);
+ }
+ return nH;
+}
+
+void
+AddToAll(THStack* all, const TH1* h, Bool_t singleStep)
+{
+ TH1* copy = static_cast<TH1*>(h->Clone(Form("%s_copy", h->GetName())));
+ copy->SetDirectory(0);
+ if (singleStep) {
+ copy->SetMarkerColor(kGray);
+ copy->SetLineColor(kGray);
+ }
+ all->Add(copy);
+}
+
+void
+DimEntry(Int_t thisId, Int_t step, TLegendEntry* e)
+{
+
+ Int_t col = (thisId == step || step <= 0) ? kBlack : kGray;
+ e->SetMarkerColor(col);
+ e->SetLineColor(col);
+ e->SetTextColor(col);
+}
+
+void
+DrawStep(THStack* deltas, THStack* nchs, THStack* prims,
+ TH1* dndeta, Int_t step)
+{
+ THStack* all = new THStack("all", "Analysis steps");
+ if (step > 0) all->SetTitle(Form("Step %d", step));
+
+ if (deltas) {
+ deltas->SetTitle("#sum_{} #Delta/#Delta_{mip}");
+ TIter next(deltas->GetHists());
+ TH1* h = 0;
+ while ((h = static_cast<TH1*>(next()))) {
+ h->SetMarkerStyle(25);
+ // Info("DrawStep", "Adding %s", h->GetName());
+ AddToAll(all, h, step>0);
+ }
+ }
+ if (nchs) {
+ nchs->SetTitle("#sum_{} N_{ch,incl}");
+ TIter next(nchs->GetHists());
+ TH1* h = 0;
+ while ((h = static_cast<TH1*>(next()))) {
+ h->SetMarkerStyle(21);
+ // Info("DrawStep", "Adding %s", h->GetName());
+ AddToAll(all, h, step>0);
+ }
+ }
+ if (prims) {
+ prims->SetTitle("#sum_{} N_{ch,primary}");
+ TIter next(prims->GetHists());
+ TH1* h = 0;
+ while ((h = static_cast<TH1*>(next()))) {
+ h->SetMarkerStyle(22);
+ // Info("DrawStep", "Adding %s", h->GetName());
+ AddToAll(all, h, step>0);
+ }
+ }
+ if (dndeta) {
+ dndeta->SetTitle("1/N dN_{ch}/d#eta");
+ dndeta->SetMarkerStyle(20);
+ dndeta->SetMarkerColor(kBlack);
+ // Info("DrawStep", "Adding %s", dndeta->GetName());
+ AddToAll(all, dndeta, step>0);
+ }
+
+ all->Draw("nostack");
+ all->GetHistogram()->SetXTitle("#eta");
+ all->GetHistogram()->SetYTitle("signal");
+ all->GetHistogram()->GetXaxis()->SetLabelFont(132);
+ all->GetHistogram()->GetXaxis()->SetTitleFont(132);
+ all->GetHistogram()->GetYaxis()->SetLabelFont(132);
+ all->GetHistogram()->GetYaxis()->SetTitleFont(132);
+ c->SetGridx();
+
+ TLegend* l = new TLegend(.33, .2, .53, .9);
+ TLegendEntry* e = 0;
+ l->SetFillColor(0);
+ l->SetFillStyle(0);
+ l->SetBorderSize(0);
+ l->SetNColumns(1);
+ l->SetTextFont(132);
+ Int_t i = 0;
+ if (deltas) {
+ TIter next(deltas->GetHists());
+ TH1* h = 0;
+ while ((h = static_cast<TH1*>(next()))) {
+ e = l->AddEntry(Form("dummy%02d", i++),h->GetTitle(),"pl");
+ e->SetMarkerStyle(20);
+ e->SetMarkerColor(h->GetMarkerColor());
+ }
+ e = l->AddEntry(Form("dummy%02d", i++), deltas->GetTitle(),"pl");
+ TH1* h = static_cast<TH1*>(deltas->GetHists()->At(0));
+ e->SetMarkerStyle(h->GetMarkerStyle());
+ DimEntry(1, step, e);
+ }
+ if (nchs) {
+ e = l->AddEntry(Form("dummy%02d",i++),nchs->GetTitle(),"pl");
+ TH1* h = static_cast<TH1*>(nchs->GetHists()->At(0));
+ e->SetMarkerStyle(h->GetMarkerStyle());
+ DimEntry(2, step, e);
+ }
+ if (prims) {
+ e = l->AddEntry(Form("dummy%02d", i++), prims->GetTitle(),"pl");
+ TH1* h = static_cast<TH1*>(prims->GetHists()->At(0));
+ e->SetMarkerStyle(h->GetMarkerStyle());
+ DimEntry(3, step, e);
+ }
+ if (dndeta) {
+ e = l->AddEntry(Form("dummy%02d", i++), dndeta->GetTitle(),"pl");
+ e->SetMarkerStyle(dndeta->GetMarkerStyle());
+ DimEntry(4, step, e);
+ }
+ l->Draw();
+
+ TString what;
+ if (step > 0) {
+ switch (step) {
+ case 1:
+ deltas->Draw("same nostack");
+ what = "After merging";
+ break;
+ case 2:
+ nchs->Draw("same nostack");
+ what = "After particle counting";
+ break;
+ case 3:
+ prims->Draw("same nostack");
+ what = "After corrections";
+ break;
+ case 4:
+ dndeta->Draw("same");
+ what = "After normalisation";
+ break;
+ default:
+ Error("DrawSteps", "Unknown step: %d (must be in 1-4)");
+ break;
+ }
+ }
+ TLatex* ltx = new TLatex(.95, .85, what);
+ ltx->SetNDC();
+ ltx->SetTextSize(.07);
+ ltx->SetTextAlign(33);
+ ltx->SetTextFont(132);
+ ltx->Draw();
+}
+
+
+void DrawSteps(const char* filename, Bool_t single)
+{
+ gStyle->SetPalette(1);
+ gStyle->SetOptFit(0);
+ gStyle->SetOptStat(0);
+
+ TFile* file = TFile::Open(filename, "READ");
+ if (!file) {
+ Error("DrawMCResult", "failed to open %s", filename);
+ return;
+ }
+ const char* fname2 = "forward_dndeta.root";
+ TFile* file2 = TFile::Open(fname2, "READ");
+ if (!file2) {
+ Error("DrawSteps", "File %s not found", fname2);
+ }
+
+ TList* forward = static_cast<TList*>(file->Get("Forward"));
+ if (!forward) {
+ Error("DrawMCResult", "List Forward not found in %s", filename);
+ return;
+ }
+ TList* forwardRes = (file2 ?
+ static_cast<TList*>(file2->Get("ForwardResults")) :
+ 0);
+ TList* forwardAll = (forwardRes ?
+ static_cast<TList*>(forwardRes->FindObject("all")) :
+ 0);
+
+
+ // THStack* res = GetStack(*forward, "ringResults", "all");
+ // THStack* mcRes = GetStack(*forward, "mcRingResults", "all");
+ THStack* deltas = GetStack(*forward, "fmdSharingFilter", "sums");
+ THStack* nchs = GetStack(*forward, "fmdDensityCalculator", "sums");
+ THStack* prims = GetStack(*forward, "fmdCorrector", "sums");
+ TH1* dndeta = (forwardAll ?
+ static_cast<TH1*>(forwardAll->FindObject("dndetaForward")):
+ 0);
+
+ Info("DrawSteps", "Got steps deltas=%p, nchs=%p, prims=%p, dndeta=%p",
+ deltas, nchs, prims, dndeta);
+
+
+ gStyle->SetTitleBorderSize(0);
+ gStyle->SetTitleFillColor(0);
+ gStyle->SetTitleStyle(0);
+ gStyle->SetTitleX(.7);
+ gStyle->SetTitleY(.95);
+ gStyle->SetTitleH(.1);
+ gStyle->SetTitleW(.25);
+ // gStyle->SetTitleColor(kBlack);
+
+
+
+ if (!single) {
+ TCanvas* c = new TCanvas("c", "C", 900, 700);
+ c->SetFillColor(0);
+ c->SetBorderSize(0);
+ c->SetTopMargin(0.05);
+ c->SetRightMargin(0.05);
+
+ DrawStep(deltas, nchs, prims, dndeta, 0);
+ return;
+ }
+ Int_t nSteps = 0;
+ if (deltas) nSteps++;
+ if (nchs) nSteps++;
+ if (prims) nSteps++;
+ if (dndeta) nSteps++;
+
+ Int_t w = (nSteps >= 4 ? 1100 : 700);
+ Int_t h = (nSteps >= 4 ? 800 : 1100);
+
+ TCanvas* c = new TCanvas("c", "C", w, h);
+ c->SetFillColor(0);
+ c->SetBorderSize(0);
+ c->SetTopMargin(0.05);
+ c->SetRightMargin(0.05);
+
+ if (nSteps >= 4)
+ c->Divide(2,(nSteps+1)/2,0,0);
+ else
+ c->Divide(1,nSteps,0,0);
+
+ for (Int_t i=1; i<=nSteps; i++) {
+ TVirtualPad* p = c->cd(i);
+ p->SetFillColor(0);
+ p->SetFillStyle(0);
+ p->SetBorderSize(0);
+ p->SetGridx();
+ p->SetGridy();
+
+ DrawStep(deltas, nchs, prims, dndeta, i);
+ }
+ c->SaveAs("steps.png");
+}
+
+
--- /dev/null
+//____________________________________________________________________
+/**
+ * Get an object with specified name from TCollection @a l
+ *
+ * @param l Collection
+ * @param name Name of object to retrieve
+ *
+ * @return Object, or null
+ */
+TObject*
+GetObject(const TObject* l, const char* name)
+{
+ if (!l->IsA()->InheritsFrom(TCollection::Class())) {
+ Error("GetObject", "passed parent %s not a TCollection, but %s",
+ l->GetName(), l->IsA()->GetName());
+ return 0;
+ }
+ TObject* o = static_cast<const TCollection*>(l)->FindObject(name);
+ if (!o) {
+ Error("GetObject", "No object '%s' found in '%s'", name, l->GetName());
+ return 0;
+ }
+ return o;
+}
+
+//____________________________________________________________________
+/**
+ * Get histogram from @a directory/which/sub/hname
+ *
+ * @param dir Directory
+ * @param which Name of parent list
+ * @param sub Name of sub-list
+ * @param hname Name of histogram
+ *
+ * @return Pointer to histogram or null
+ */
+TH1D*
+GetHist(TDirectory* dir,
+ const char* which,
+ const char* sub,
+ const char* hname)
+{
+ if (!dir) return 0;
+ TList* parent = static_cast<TList*>(dir->Get(which));
+ if (!parent) {
+ Error("GetHist", "List '%s' not found in '%s'", which, dir->GetName());
+ return 0;
+ }
+ TList* child = static_cast<TList*>(parent->FindObject(sub));
+ if (!child) {
+ Error("GetHist", "List '%s' not found in '%s'", sub, parent->GetName());
+ return 0;
+ }
+ TObject* o = GetObject(child,hname);
+ if (!o) return 0;
+ return static_cast<TH1D*>(o);
+}
+
+
+//____________________________________________________________________
+/**
+ * Get histogram from
+ * <i>dir/which</i>Results<i>/sub/</i>dndeta<i>which</i>[_rebin<i>rebin</i>]
+ *
+ * @param dir Directory
+ * @param which Name
+ * @param rebin Optional rebinning
+ * @param sub Sub-list name
+ *
+ * @return Histogram pointer or null
+ */
+TH1D*
+GetHist(TDirectory* dir,
+ const char* which,
+ UShort_t rebin,
+ const char* sub="all")
+{
+ TString name(Form("dndeta%s", which));
+ if (rebin > 1) name.Append(Form("_rebin%02d", rebin));
+ return GetHist(dir, Form("%sResults", which), sub, name);
+}
+
+//____________________________________________________________________
+/**
+ * Merge two histograms into one
+ *
+ * @param cen Central part
+ * @param fwd Forward part
+ * @param xlow On return, lower eta bound
+ * @param xhigh On return, upper eta bound
+ *
+ * @return Newly allocated histogram or null
+ */
+TH1*
+Merge(const TH1* cen, const TH1* fwd, Double_t& xlow, Double_t& xhigh)
+{
+ TH1* tmp = static_cast<TH1*>(fwd->Clone("dndetaMerged"));
+ // tmp->SetMarkerStyle(28);
+ // tmp->SetMarkerColor(kBlack);
+ tmp->SetDirectory(0);
+ xlow = 100;
+ xhigh = -100;
+ for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
+ Double_t cc = cen->GetBinContent(i);
+ Double_t cf = fwd->GetBinContent(i);
+ Double_t ec = cen->GetBinError(i);
+ Double_t ef = fwd->GetBinError(i);
+ Double_t nc = cf;
+ Double_t ne = ef;
+ if (cc < 0.001 && cf < 0.01) continue;
+ xlow = TMath::Min(tmp->GetXaxis()->GetBinLowEdge(i),xlow);
+ xhigh = TMath::Max(tmp->GetXaxis()->GetBinUpEdge(i),xhigh);
+ if (cc > 0.001) {
+ nc = cc;
+ ne = ec;
+ if (cf > 0.001) {
+ nc = (cf + cc) / 2;
+ ne = TMath::Sqrt(ec*ec + ef*ef);
+ }
+ }
+ tmp->SetBinContent(i, nc);
+ tmp->SetBinError(i, ne);
+ }
+ return tmp;
+}
+
+//____________________________________________________________________
+/**
+ * Function to calculate
+ * @f[
+ * g(x;A_1,A_2,\sigma_1,\sigma_2) =
+ * A_1\left(\frac{1}{2\pi\sigma_1}e^{(x/\sigma_1)^2} -
+ * A_2\frac{1}{2\pi\sigma_2}e^{(x/\sigma_2)^2}\right)
+ * @f]
+ *
+ * @param xp Pointer to x array
+ * @param pp Pointer to parameter array
+ *
+ * @return @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$
+ */
+Double_t myFunc(Double_t* xp, Double_t* pp)
+{
+ Double_t x = xp[0];
+ Double_t a1 = pp[0];
+ Double_t a2 = pp[1];
+ Double_t s1 = pp[2];
+ Double_t s2 = pp[3];
+ return a1*(TMath::Gaus(x, 0, s1) - a2 * TMath::Gaus(x, 0, s2));
+}
+//____________________________________________________________________
+/**
+ * Calculate
+ * @f[
+ * r(x) = \frac{g(x;A_1,A_2,\sigma_1,\sigma_2)}{
+ * g(x;A_1',A_2',\sigma'_1,\sigma'_2)}
+ * @f]
+ *
+ * @param xp Pointer to X array
+ * @param pp Pointer to parameter array (8 entries)
+ *
+ * @return @f$r(x)@f$
+ */
+Double_t myRatio(Double_t* xp, Double_t* pp)
+{
+ Double_t denom = myFunc(xp, &(pp[4]));
+ if (TMath::Abs(denom) < 1.e-6) return 0;
+ return myFunc(xp, pp) / denom;
+}
+
+//____________________________________________________________________
+/**
+ * Fit @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$ to histogram data
+ *
+ * @param tmp Histogram
+ * @param xlow Lower x bound
+ * @param xhigh Upper x bound
+ *
+ * @return Fitted function
+ */
+TF1*
+FitMerged(TH1* tmp, Double_t xlow, Double_t xhigh)
+{
+ TF1* tmpf = new TF1("tmpf", "gaus", xlow, xhigh);
+ tmp->Fit(tmpf, "N", "");
+ tmp->SetDirectory(0);
+
+ TF1* fit = new TF1("f", myFunc, xlow, xhigh, 4);
+ fit->SetParNames("a_{1}", "a_{2}", "#sigma_{1}", "#sigma_{2}");
+ fit->SetParameters(tmpf->GetParameter(0),
+ .2,
+ tmpf->GetParameter(2),
+ tmpf->GetParameter(2)/4);
+ fit->SetParLimits(3, 0, 100);
+ fit->SetParLimits(4, 0, 100);
+ tmp->Fit(fit,"0W","");
+ return fit;
+}
+
+//____________________________________________________________________
+/**
+ * Make band of systematic errors
+ *
+ * @param tmp Histogram
+ * @param fit Fit
+ */
+void
+MakeSysError(TH1* tmp, TF1* fit)
+{
+ for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
+ Double_t tc = tmp->GetBinContent(i);
+ if (tc < 0.01) continue;
+ Double_t x = tmp->GetXaxis()->GetBinCenter(i);
+ Double_t y = fit->Eval(x);
+ tmp->SetBinContent(i, y);
+ tmp->SetBinError(i,sysErr*y);
+ }
+ return tmp;
+}
+
+//____________________________________________________________________
+/**
+ * Transform a graph into a histogram
+ *
+ * @param g
+ *
+ * @return
+ */
+TH1*
+Graph2Hist(const TGraphAsymmErrors* g)
+{
+ Int_t nBins = g->GetN();
+ TArrayF bins(nBins+1);
+ TArrayF y(nBins);
+ TArrayF ey(nBins);
+ Double_t dx = 0;
+ Double_t xmin = 10000;
+ Double_t xmax = -10000;
+ for (Int_t i = 0; i < nBins; i++) {
+ Double_t x = g->GetX()[i];
+ Double_t exl = g->GetEXlow()[i];
+ Double_t exh = g->GetEXhigh()[i];
+ xmin = TMath::Min(x-exl, xmin);
+ xmax = TMath::Max(x+exh, xmax);
+ bins.fArray[i] = x-exl;
+ bins.fArray[i+1] = x+exh;
+ Double_t dxi = exh+exl;
+ if (dxi == 0 && i != 0) dxi = bins.fArray[i]-bins.fArray[i-1];
+ if (dx == 0) dx = dxi;
+ else if (dxi != dx) dx = 0;
+
+ y.fArray[i] = g->GetY()[i];
+ ey.fArray[i] = TMath::Max(g->GetEYlow()[i],g->GetEYhigh()[i]);
+
+ }
+ TString name(g->GetName());
+ TString title(g->GetTitle());
+ TH1D* h = 0;
+ if (dx != 0) {
+ h = new TH1D(name.Data(), title.Data(), nBins,
+ bins[0]-dx/2, bins[nBins]+dx/2);
+ }
+ else {
+ h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray);
+ }
+ for (Int_t i = 1; i <= nBins; i++) {
+ h->SetBinContent(i, y.fArray[i-1]);
+ h->SetBinError(i, ey.fArray[i-1]);
+ }
+ h->SetMarkerStyle(g->GetMarkerStyle());
+ h->SetMarkerColor(g->GetMarkerColor());
+ h->SetMarkerSize(g->GetMarkerSize());
+ h->SetDirectory(0);
+
+ return h;
+}
+
+//____________________________________________________________________
+/**
+ * Calculate ratio of histogram to function
+ *
+ * @param h Histogram
+ * @param f Function
+ * @param title (Optional) title
+ *
+ * @return Ratio in a histogram
+ */
+TH1*
+Ratio(TH1* h, TF1* f, const char* title)
+{
+ TH1* ret = static_cast<TH1*>(h->Clone(Form("%s_%s",
+ h->GetName(),
+ f->GetName())));
+ ret->SetDirectory(0);
+ if (title) ret->SetTitle(title);
+ else ret->SetTitle(Form("%s (data) / %s",
+ h->GetTitle(),f->GetTitle()));
+ ret->Reset();
+ for (Int_t i = 1; i <= ret->GetNbinsX(); i++) {
+ Double_t cc = h->GetBinContent(i);
+ if (cc < 0.01) {
+ ret->SetBinContent(i,0);
+ ret->SetBinError(i,0);
+ continue;
+ }
+ Double_t xx = h->GetBinCenter(i);
+ Double_t ee = h->GetBinError(i);
+ Double_t ff = f->Eval(xx);
+ Double_t yy = cc / ff;
+ Double_t ey = ee / ff;
+ ret->SetBinContent(i, yy);
+ ret->SetBinError(i, ey);
+ }
+ return ret;
+}
+
+//____________________________________________________________________
+/**
+ * Get the UA5 data
+ *
+ * @param type Trigger type (1: INEL, 4: NSD)
+ * @param p On return, positive part or null
+ * @param n On return, negative part or null
+ * @param xlow On return, lower X bound
+ * @param xhigh On return, upper X bound
+ *
+ * @return Merged histogram or null
+ */
+TH1D*
+GetUA5Data(UShort_t type, TH1*& p, TH1*& n,
+ Double_t& xlow, Double_t& xhigh)
+{
+ gROOT->SetMacroPath(Form(".:$ALICE_ROOT.trunk/PWG2/FORWARD/analysis2/:%s",
+ gROOT->GetMacroPath()));
+ gROOT->LoadMacro("OtherData.C");
+
+ p = 0;
+ n = 0;
+ TGraphAsymmErrors* gp = GetSingle(UA5, 1, 900, type, 0, 0);
+ TGraphAsymmErrors* gn = GetSingle(UA5+10, 1, 900, type, 0, 0);
+ if (!gp || !gn) return 0;
+
+ p = Graph2Hist(gp);
+ n = Graph2Hist(gn);
+
+ Int_t nn = p->GetNbinsX();
+ xlow = n->GetXaxis()->GetXmin();
+ xhigh = p->GetXaxis()->GetXmax();
+ TH1D* ret = new TH1D("ua5", "UA5", 2*nn, xlow, xhigh);
+ ret->SetDirectory(0);
+ ret->SetMarkerColor(p->GetMarkerColor());
+ ret->SetMarkerStyle(p->GetMarkerStyle());
+
+ for (Int_t i = 1; i <= nn; i++) {
+ ret->SetBinContent(nn+i, p->GetBinContent(i));
+ ret->SetBinContent( i, n->GetBinContent(i));
+ ret->SetBinError(nn+i, p->GetBinError(i));
+ ret->SetBinError( i, n->GetBinError(i));
+ }
+ return ret;
+}
+
+
+//____________________________________________________________________
+/**
+ *
+ *
+ */
+void
+DrawUA5Ratios(const char* fname="forward_dndeta.root", UShort_t rebin=5)
+{
+ TFile* forward_dndeta = TFile::Open(fname, "READ");
+ if (!forward_dndeta) {
+ Error("DrawUA5Ratios", "%s not found", fname);
+ return;
+ }
+
+ UShort_t type = 1;
+ TH1D* forward = GetHist(forward_dndeta, "Forward", rebin);
+ TH1D* central = GetHist(forward_dndeta, "Central", rebin);
+
+ TObject* sys = GetObject(forward_dndeta->Get("ForwardResults"), "sys");
+ TObject* sNN = GetObject(forward_dndeta->Get("ForwardResults"), "sNN");
+ TObject* trg = GetObject(forward_dndeta->Get("ForwardResults"), "trigger");
+ if (!sys || !sNN || !trg) return;
+ if (sys->GetUniqueID() != 1) {
+ Error("DrawUA5Ratios", "Comparison only valid for pp, not %s",
+ sys->GetTitle());
+ return;
+ }
+ if (sNN->GetUniqueID() != 900) {
+ Error("DrawUA5Ratios", "Comparison only valid for 900GeV, not %dGeV",
+ sNN->GetUniqueID());
+ return;
+ }
+ if (trg->GetUniqueID() != 1 &&
+ trg->GetUniqueID() != 4) {
+ Error("DrawUA5Ratios",
+ "Comparison only valid for INEL or NSD, not %s (%d)",
+ trg->GetTitle(), trg->GetUniqueID());
+ return;
+ }
+
+
+ Double_t alilow, alihigh;
+ TH1D* ali = Merge(central, forward, alilow, alihigh);
+ TF1* alifit = FitMerged(ali, alilow, alihigh);
+ ali->SetTitle("Forward+Central");
+ alifit->SetLineColor(kRed+1);
+ alifit->SetLineStyle(2);
+ alifit->SetName("alifit");
+ alifit->SetTitle("Fit to Forward+Central");
+
+ Double_t ua5low, ua5high;
+ TH1* ua5p, *ua5n;
+ TH1D* ua5 = GetUA5Data(trg->GetUniqueID(), ua5p, ua5n, ua5low, ua5high);
+ TF1* ua5fit = FitMerged(ua5, ua5low, ua5high);
+ ua5fit->SetLineColor(kBlue+1);
+ ua5fit->SetLineStyle(3);
+ ua5fit->SetName("ua5fit");
+ ua5fit->SetTitle("Fit to UA5");
+
+ gStyle->SetOptTitle(0);
+ TCanvas* c = new TCanvas("c", "C", 900, 900);
+ c->SetFillColor(0);
+ c->SetFillStyle(0);
+ c->SetBorderMode(0);
+ c->SetBorderSize(0);
+
+ Double_t yd = .3;
+
+ TPad* p1 = new TPad("p1", "p1", 0, yd, 1, 1, 0, 0, 0);
+ p1->SetBorderSize(0);
+ p1->SetBorderMode(0);
+ p1->SetFillColor(0);
+ p1->SetFillStyle(0);
+ p1->SetRightMargin(0.02);
+ p1->SetTopMargin(0.02);
+ p1->SetBottomMargin(0.00);
+ p1->SetGridx();
+ p1->Draw();
+ p1->cd();
+
+ THStack* all = new THStack("all", "All");
+ all->Add(ua5p);
+ all->Add(ua5n);
+ // all->Add(ua5);
+ all->Add(forward);
+ all->Add(central);
+ // all->Add(merged);
+ all->Draw("nostack");
+ all->SetMinimum(-.07);
+ all->GetXaxis()->SetTitleFont(132);
+ all->GetYaxis()->SetTitleFont(132);
+ all->GetXaxis()->SetLabelFont(132);
+ all->GetYaxis()->SetLabelFont(132);
+ all->GetYaxis()->SetDecimals();
+ p1->Clear();
+ all->Draw("nostack");
+ // ua5p->Draw("same p");
+ // ua5m->Draw("same p");
+ alifit->Draw("same");
+ ua5fit->Draw("same");
+
+ TLegend* l = new TLegend(.2, .1, .8, .5,
+ Form("pp @ #sqrt{s}=900GeV, %s",trg->GetTitle()));
+ l->AddEntry(ua5, Form("U: %s", ua5->GetTitle()), "lp");
+ l->AddEntry(forward, "F: Forward", "lp");
+ l->AddEntry(central, "C: Central", "lp");
+ l->AddEntry(alifit,
+ Form("f: %s: %4.2f#left[e^{(x/%4.2f)^{2}}-"
+ "%4.2fe^{(x/%4.2f)^{2}}#right]",
+ alifit->GetTitle(),
+ alifit->GetParameter(0),
+ alifit->GetParameter(2),
+ alifit->GetParameter(1),
+ alifit->GetParameter(3)), "l");
+ l->AddEntry(ua5fit,
+ Form("u: %s: %4.2f#left[e^{(x/%4.2f)^{2}}-"
+ "%4.2fe^{(x/%4.2f)^{2}}#right]",
+ ua5fit->GetTitle(),
+ ua5fit->GetParameter(0),
+ ua5fit->GetParameter(2),
+ ua5fit->GetParameter(1),
+ ua5fit->GetParameter(3)), "l");
+ l->SetTextFont(132);
+ l->SetBorderSize(0);
+ l->SetFillColor(0);
+ l->SetFillStyle(0);
+ l->Draw();
+
+ c->cd();
+ TPad* p2 = new TPad("p2", "p2", 0, 0, 1, yd, 0, 0, 0);
+ p2->SetBorderSize(0);
+ p2->SetBorderMode(0);
+ p2->SetFillColor(0);
+ p2->SetFillStyle(0);
+ p2->SetRightMargin(0.02);
+ p2->SetTopMargin(0.00);
+ p2->SetBottomMargin(0.15);
+ p2->SetGridx();
+ p2->Draw();
+ p2->cd();
+
+ THStack* ratios = new THStack("Ratios", "Ratios");
+ TH1* ua5ali = Ratio(ua5, alifit, 0);
+ TH1* aliua5 = Ratio(ali, ua5fit, 0);
+ ratios->Add(ua5ali);
+ ratios->Add(aliua5);
+ ratios->Draw("nostack");
+ ratios->SetMinimum(0.4);
+ ratios->GetYaxis()->SetTitleSize(2*ratios->GetYaxis()->GetTitleSize());
+ ratios->GetYaxis()->SetLabelSize(2*ratios->GetYaxis()->GetLabelSize());
+ ratios->GetYaxis()->SetNdivisions(508);
+ ratios->GetXaxis()->SetTitleSize(2*ratios->GetXaxis()->GetTitleSize());
+ ratios->GetXaxis()->SetLabelSize(2*ratios->GetXaxis()->GetLabelSize());
+ ratios->GetXaxis()->SetNdivisions(510);
+ ratios->GetXaxis()->SetTitle("#eta");
+ ratios->GetXaxis()->SetTitleFont(132);
+ ratios->GetYaxis()->SetTitleFont(132);
+ ratios->GetXaxis()->SetLabelFont(132);
+ ratios->GetYaxis()->SetLabelFont(132);
+ ratios->GetYaxis()->SetDecimals();
+ p2->Clear();
+
+ TGraphErrors* sysErr = new TGraphErrors(2);
+ sysErr->SetPoint(0, all->GetHistogram()->GetXaxis()->GetXmin(),1);
+ sysErr->SetPoint(1, all->GetHistogram()->GetXaxis()->GetXmax(),1);
+ sysErr->SetPointError(0,0,.07);
+ sysErr->SetPointError(1,0,.07);
+ sysErr->SetTitle("Systematic error on Forward+Central");
+ sysErr->SetFillColor(kYellow+1);
+ sysErr->SetFillStyle(3001);
+ sysErr->SetHistogram(ratios->GetHistogram());
+ sysErr->DrawClone("ael3");
+ ratios->DrawClone("nostack same");
+
+ TF1* fitfit = new TF1("fitfit", myRatio, alilow, alihigh, 8);
+ fitfit->SetParameters(ua5fit->GetParameter(0),
+ ua5fit->GetParameter(1),
+ ua5fit->GetParameter(2),
+ ua5fit->GetParameter(3),
+ alifit->GetParameter(0),
+ alifit->GetParameter(1),
+ alifit->GetParameter(2),
+ alifit->GetParameter(3));
+ fitfit->Draw("same");
+
+ TLegend* ll = new TLegend(.3,.15, .7, .6);
+ ll->AddEntry(sysErr, sysErr->GetTitle(), "f");
+ ll->AddEntry(ua5ali, ua5ali->GetTitle(), "lp");
+ ll->AddEntry(aliua5, aliua5->GetTitle(), "lp");
+ ll->AddEntry(fitfit, "UA5 (fit) / Forward+Central (fit)", "lp");
+ ll->SetTextFont(132);
+ ll->SetBorderSize(0);
+ ll->SetFillColor(0);
+ ll->SetFillStyle(0);
+ ll->Draw();
+
+
+ c->SaveAs(Form("ua5_ratios_%s_r%02d.png", trg->GetTitle(), rebin));
+ gROOT->GetInterpreter()->UnloadFile("OtherData.C");
+
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
# include <THStack.h>
# include <TLegend.h>
# include <TLegendEntry.h>
+# include <TLatex.h>
# include <TLine.h>
# include <TString.h>
# include <TCanvas.h>
# include <TError.h>
+# include <TColor.h>
#else
class THStack;
class TAxis;
return static_cast<THStack*>(o);
}
+TH1*
+GetHist(const TList* list, const char* name, Int_t rebin)
+{
+ if (!list) {
+ Warning("GetStack", "List is null");
+ return 0;
+ }
+ TList* all = static_cast<TList*>(list->FindObject("all"));
+ if (!all) {
+ Warning("GetHist", "List all not found in %s", list->GetName());
+ return 0;
+ }
+
+ TString n(name);
+ if (rebin > 1) n.Append(Form("_rebin%02d", rebin));
+ TObject* o = all->FindObject(n);
+ if (!o) {
+ Warning("GetStack", "No %s object found in %s", n.Data(), all->GetName());
+ return 0;
+ }
+ return static_cast<TH1*>(o);
+}
+
/**
* Add histograms from one stack to another
*
void
BuildLegend(const THStack* stack, const TAxis* c)
{
- TLegend* l = new TLegend(.75, .8, 1-gPad->GetRightMargin(),
- 1-gPad->GetTopMargin(), "");
+ Double_t x1 = .75, x2 = 1-gPad->GetRightMargin();
+ Double_t y1 = .8, y2 = 1-gPad->GetTopMargin();
+ if (!c) {
+ // PP
+ x1 = .4;
+ y1 = .4;
+ x2 = .8;
+ y2 = .6;
+ }
+ TLegend* l = new TLegend(x1, y1, x2, y2, "");
l->SetFillColor(0);
l->SetFillStyle(0);
l->SetBorderSize(0);
if (c) dd->SetMarkerColor(kBlack);
else dd->SetMarkerColor(color);
dd->SetMarkerStyle(style);
+ Int_t st = dd->GetMarkerStyle();
+ if (st == 27 || st == 28 || st == 29 || st == 30 || st == 33 || st == 34)
+ dd->SetMarkerSize(1.4*dd->GetMarkerSize());
}
// TLegendEntry* sep = l->AddEntry("s", "_", "h");
// sep->SetTextSize(0.01);
TLine* sep = new TLine(0,0,1,1);
sep->SetLineWidth(1);
- sep->DrawLineNDC(.77, .795, 1-gPad->GetRightMargin()-.03, .795);
+ sep->DrawLineNDC(x1+.02, y1-.005, x2-.03, y1-.005);
// sep->Draw();
- TLegend* l2 = new TLegend(.75, .7, 1-gPad->GetRightMargin(), .79, "");
+ TLegend* l2 = new TLegend(x1, y1-.005, x2, y1-.15, "");
l2->SetFillColor(0);
l2->SetFillStyle(0);
l2->SetBorderSize(0);
d2->SetLineColor(kBlack);
d2->SetMarkerColor(kBlack);
d2->SetMarkerStyle(24);
-
+
l->Draw();
l2->Draw();
}
+
+void
+AddInformation(TList* forward, bool prelim=true)
+{
+ Double_t x = .12;
+ Double_t y = .95;
+ Double_t ts = 0.05;
+ if (prelim) {
+ TLatex* wip = new TLatex(x, y, "Work in progress");
+ wip->SetNDC();
+ wip->SetTextColor(TColor::GetColor(234,26,46));
+ wip->SetTextAlign(13);
+ wip->SetTextFont(132);
+ wip->SetTextSize(ts);
+ wip->Draw();
+ y -= ts;
+ }
+
+ TObject* sNN = forward->FindObject("sNN");
+ TObject* sys = forward->FindObject("sys");
+ TObject* trg = forward->FindObject("trigger");
+ TObject* vtx = forward->FindObject("vtxAxis");
+ TObject* sch = forward->FindObject("scheme");
+
+ TString t(sys->GetTitle());
+ Bool_t isPP = t == "pp";
+
+ TString s = Form("%s @ #sqrt{s%s}=",
+ sys->GetTitle(),
+ (isPP ? "" : "_{NN}"));
+ Int_t cms = sNN->GetUniqueID();
+ if (cms > 1000) s.Append(Form("%5.2fTeV", float(cms)/1000));
+ else s.Append(Form("%3dGeV", cms));
+ s.Append(Form(", %s", trg->GetTitle()));
+
+ // if (isPP) { x = .3; y = .3; }
+ if (isPP) {
+ x = 1-gPad->GetRightMargin()-.01;
+ y = 1-gPad->GetTopMargin()-.01;
+ ts = .04;
+ }
+ TLatex* ltx = new TLatex(x, y, s.Data());
+ ltx->SetNDC();
+ ltx->SetTextColor(TColor::GetColor(41,73,156));
+ ltx->SetTextAlign((isPP ? 33 : 13));
+ ltx->SetTextFont(132);
+ ltx->SetTextSize(ts);
+ ltx->Draw();
+ y -= ts;
+ ltx->DrawLatex(x, y, vtx->GetTitle());
+ y -= ts;
+ ltx->DrawLatex(x, y, sch->GetTitle());
+}
+
+Double_t myFunc(Double_t* xp, Double_t* pp)
+{
+ Double_t x = xp[0];
+ Double_t a1 = pp[0];
+ Double_t a2 = pp[1];
+ Double_t s1 = pp[2];
+ Double_t s2 = pp[3];
+ return a1*(TMath::Gaus(x, 0, s1) - a2 * TMath::Gaus(x, 0, s2));
+}
+
+TH1*
+MakeSysError(const TH1* cen, const TH1* fwd, Double_t sysErr=0.7)
+{
+ TH1* tmp = static_cast<TH1*>(fwd->Clone("dndetaFitted"));
+ tmp->SetMarkerStyle(0);
+ tmp->SetFillColor(kGray);
+ tmp->SetFillStyle(3001);
+ tmp->SetDirectory(0);
+ Double_t xlow = 100;
+ Double_t xhigh = 0;
+ for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
+ Double_t cc = cen->GetBinContent(i);
+ Double_t cf = fwd->GetBinContent(i);
+ Double_t ec = fwd->GetBinError(i);
+ Double_t ef = fwd->GetBinError(i);
+ Double_t nc = cf;
+ Double_t ne = ef;
+ if (cc < 0.001 && cf < 0.01) continue;
+ xlow = TMath::Min(tmp->GetXaxis()->GetBinLowEdge(i),xlow);
+ xhigh = TMath::Max(tmp->GetXaxis()->GetBinUpEdge(i),xhigh);
+ if (cc > 0.001) {
+ nc = cc;
+ ne = ec;
+ if (cf > 0.001) {
+ nc = (cf + cc) / 2;
+ ne = TMath::Sqrt(ec*ec + ef*ef);
+ }
+ }
+ tmp->SetBinContent(i, nc);
+ tmp->SetBinError(i, ne);
+ }
+ TF1* tmpf = new TF1("tmpf", "gaus", xlow, xhigh);
+ tmp->Fit(tmpf, "NQ", "");
+
+ TF1* fit = new TF1("f", myFunc, xlow, xhigh, 4);
+ fit->SetParNames("a_{1}", "a_{2}", "#sigma_{1}", "#sigma_{2}");
+ fit->SetParameters(tmpf->GetParameter(0),
+ .2,
+ tmpf->GetParameter(2),
+ tmpf->GetParameter(2)/4);
+ tmp->Fit(fit,"0WQ","");
+ for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
+ Double_t tc = tmp->GetBinContent(i);
+ if (tc < 0.01) continue;
+ Double_t x = tmp->GetXaxis()->GetBinCenter(i);
+ Double_t y = fit->Eval(x);
+ tmp->SetBinContent(i, y);
+ tmp->SetBinError(i,sysErr*y);
+ }
+ return tmp;
+}
/**
* Function to draw the results from forward_dndeta.root file
* @param filename File to open and draw stuff from >
*/
void
-SimpledNdeta(Int_t rebin=5, const char* filename="forward_dndeta.root")
+SimpledNdeta(Int_t what=0x5,
+ Int_t rebin=5, const char* filename="forward_dndeta.root")
{
TFile* file = TFile::Open(filename, "READ");
if (!file) {
TList* forward = static_cast<TList*>(file->Get("ForwardResults"));
TList* central = static_cast<TList*>(file->Get("CentralResults"));
+ TList* mctruth = static_cast<TList*>(file->Get("MCTruthResults"));
THStack* all = new THStack("all", "All");
- AddStack(all, forward, "dndeta", rebin);
- AddStack(all, forward, "dndetaMC", rebin);
- AddStack(all, forward, "dndetaTruth", rebin);
- AddStack(all, central, "dndeta", rebin);
+ if (what & 0x1) AddStack(all, forward, "dndeta", rebin);
+ if (what & 0x2) AddStack(all, forward, "dndetaMC", rebin);
+ if (what & 0x1) AddStack(all, central, "dndeta", rebin);
+ if (what & 0x2) AddStack(all, central, "dndetaMC", rebin);
+ if (what & 0x4) AddStack(all, mctruth, "dndeta", rebin);
+
+ TH1* tmp = 0;
+ if (what & 0x1) {
+ Double_t sysErr = 0.07;
+ TH1* fwd = GetHist(forward, "dndetaForward", rebin);
+ TH1* cen = GetHist(central, "dndetaCentral", rebin);
+ tmp = MakeSysError(cen, fwd, sysErr);
+ // fit = tmpf;
+ }
TAxis* centAxis = static_cast<TAxis*>(forward->FindObject("centAxis"));
+ Bool_t isPP = centAxis == 0;
gStyle->SetPalette(1);
gStyle->SetOptTitle(0);
ya->SetLabelFont(132);
ya->SetTitle("dN_{ch}/d#eta");
+ if (tmp) {
+ tmp->Draw("e5 same");
+ all->Draw("same nostack");
+ }
+
+
+ // if (fit) fit->Draw("same");
+ // if (tmp) tmp->Draw("same");
+
BuildCentLegend(centAxis);
BuildLegend(all, centAxis);
- c->SaveAs("dndeta.C");
+ AddInformation(forward);
+
+ c->SaveAs("dndeta_simple.C");
+ c->SaveAs("dndeta_simple.png");
+ c->SaveAs("dndeta_simple.root");
}
//
// EOF
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff