#include <TF1.h>
#include <TSpectrum.h>
#include <TLegend.h>
+#include <TLatex.h>
#include <TLine.h>
+#include <TPolyMarker.h>
+#include <TSpectrum.h>
+#include <TList.h>
+#include <TGraph.h>
+
+Double_t landau(Double_t* xp, Double_t* pp)
+{
+ Double_t x = xp[0];
+ Double_t A = pp[0];
+ Double_t mpv = pp[1];
+ Double_t w = pp[2];
+
+ Double_t v = mpv; // + w * 0.22278;
+ return A * TMath::Landau(x, v, w);
+}
+
+Double_t foldLandau(Double_t* xp, Double_t* pp)
+{
+ Double_t x = xp[0];
+ Double_t A = pp[0];
+ Double_t mpv = pp[1];
+ Double_t w = pp[2];
+ Double_t B = pp[3];
+ Double_t sigma = pp[4];
+
+ return A * (TMath::Landau(x,mpv,w) + B * TMath::Gaus(x, mpv, sigma));
+}
/** @class DrawESD
@brief Draw digit ADC versus Rec point mult
private:
TH1D* fMult; // Histogram
const Double_t fCorr;
+ TList fCleanup;
public:
//__________________________________________________________________
- DrawESD(Int_t n=1000, Double_t mmin=-0.5, Double_t mmax=20.5)
+ DrawESD(Int_t n=2000, Double_t mmin=-0.5, Double_t mmax=15.5)
: fCorr(1) // 0.68377 / 1.1)
{
AddLoad(kESD);
- fMult = new TH1D("mult", " Multiplicity (strip)", n, mmin, mmax);
+ fMult = new TH1D("mult", "#DeltaE/#DeltaE_{MIP)", n, mmin, mmax);
fMult->Sumw2();
- fMult->SetXTitle("Strip Multiplicity");
+ fMult->SetFillColor(kRed+1);
+ fMult->SetFillStyle(3001);
+ fMult->SetXTitle("#DeltaE/#DeltaE_{MIP}");
+ fCleanup.Add(fMult);
+ }
+ ~DrawESD()
+ {
+ fCleanup.Delete();
}
//__________________________________________________________________
/** Begining of event
return AliFMDInput::Begin(ev);
}
//__________________________________________________________________
- Bool_t ProcessESD(UShort_t det,
- Char_t rng,
- UShort_t sec,
- UShort_t str,
- Float_t /* eta */,
+ Bool_t ProcessESD(UShort_t /* det */,
+ Char_t /* rng */,
+ UShort_t /* sec */,
+ UShort_t /* str */,
+ Float_t eta,
Float_t mult)
{
// Cache the energy loss
- if (mult > 0)
- Info("ProcessESD", "FMD%d%c[%2d,%3d]=%f", det, rng, sec, str, mult);
- if (mult/fCorr > 0.001) fMult->Fill(mult/fCorr);
+ // if (mult > 0)
+ // Info("ProcessESD", "FMD%d%c[%2d,%3d]=%f", det, rng, sec, str, mult);
+ if (mult <= 0 || mult == AliESDFMD::kInvalidMult) return kTRUE;
+ Double_t x = mult;
+ if (!fESD->IsAngleCorrected()) {
+ Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
+ Double_t corr = TMath::Abs(TMath::Cos(theta));
+ Double_t cmult = corr * mult;
+ x = cmult;
+ }
+ if (x > 0.001) fMult->Fill(x);
+ return kTRUE;
+ }
+ //__________________________________________________________________
+ Bool_t ProcessESDs()
+ {
+ if (!AliFMDInput::ProcessESDs()) return kFALSE;
+ // Info("ProcessESDs", "ESD is %sangle corrected",
+ // fESD->IsAngleCorrected() ? "" : "not ");
return kTRUE;
}
//__________________________________________________________________
TF1* FitPeak(Int_t n, TH1D* hist, Double_t min, Double_t& max)
{
if (TMath::Abs(max-min) < .25) return 0;
- std::cout << "Fit peack in range " << min << " to " << max << std::endl;
+ std::cout << "Fit peak in range " << min << " to " << max << std::endl;
TF1* l = new TF1(Form("l%02d", n), "landau", min, max);
hist->GetXaxis()->SetRangeUser(0, 4);
hist->Fit(l, "0Q", "", min, max);
var = hist->GetRMS();
}
+ //__________________________________________________________________
+ const char* PrettyFloat(float x)
+ {
+ if (x == 0) return Form("%9.4f", x);
+ float e = TMath::Floor(TMath::Log10(TMath::Abs(x)));
+ if (TMath::Abs(e) < 4) {
+ return Form("%9.4f", x);
+ }
+ float f = x / TMath::Power(10,e);
+ return Form("%4.2f#times10^{%d}", f, int(e));
+ }
+ //__________________________________________________________________
+ void ShowFit(Double_t x1, Double_t y1, const char* title,
+ TF1* f, Double_t dx=0, Double_t dy=0.05)
+ {
+ Double_t x = x1, y = y1;
+ TLatex* latex = new TLatex(x, y, title);
+ latex->SetTextFont(132);
+ latex->SetTextSize(0.8*dy);
+ latex->SetNDC();
+ latex->Draw();
+ x -= dx;
+ y -= dy;
+ const Double_t eqDx=0.1;
+ Double_t chi2 = f->GetChisquare();
+ Int_t ndf = f->GetNDF();
+ Double_t prob = f->GetProb();
+ latex->DrawLatex(x, y, "#chi^{2}/NDF");
+ latex->DrawLatex(x+eqDx, y, Form("= %7.4f/%3d=%5.2f (%7.3f%%)",
+ chi2, ndf, chi2/ndf, 100*prob));
+ Int_t n = f->GetNpar();
+ Double_t* p = f->GetParameters();
+ Double_t* e = f->GetParErrors();
+ for (int i = 0; i < n; i++) {
+ x -= dx;
+ y -= dy;
+ latex->DrawLatex(x, y, f->GetParName(i));
+ latex->DrawLatex(x+eqDx, y, Form("= %s", PrettyFloat(p[i])));
+ latex->DrawLatex(x+2.2*eqDx, y, Form("#pm %s", PrettyFloat(e[i])));
+ }
+ }
//__________________________________________________________________
Bool_t Finish()
{
if (fMult->GetEntries() <= 0) return kFALSE;
- TCanvas* c = new TCanvas("c", "C");
+ TCanvas* c = new TCanvas("c", "C", 1200, 800);
+ fCleanup.Add(c);
c->cd();
- // c->SetLogy();
- fMult->GetXaxis()->SetRangeUser(0,4);
- fMult->Scale(1. / fMult->GetEntries());
- fMult->SetStats(kFALSE);
- fMult->SetFillColor(2);
- fMult->SetFillStyle(3001);
- fMult->Draw("hist e");
+ c->SetLogy();
+ c->SetTopMargin(0.05);
+ c->SetRightMargin(0.05);
+ c->SetFillColor(0);
+ c->SetBorderMode(0);
+
+ TLegend* leg = new TLegend(.1, .1, .4, .2);
+ leg->SetFillColor(0);
+ leg->SetBorderSize(1);
+
+ DrawMult(c, leg);
+
+ Double_t xmax=0, xmin=0, ymax=0;
+ FindMinMax(xmin, xmax, ymax);
+
+ FitLandau(xmin, xmax, ymax, leg);
+ DrawResponse(xmax, ymax, leg);
+ // TF1* f = FitMultiLandau(xmin, xmax, leg);
+ // DrawLandaus(f);
+
+ c->cd();
+ leg->Draw();
+
+ c->Modified();
+ c->Update();
+ c->cd();
+ c->SaveAs("esd_eloss.png");
+
+ fCleanup.Add(leg);
+
return kTRUE;
- Double_t mean, rms;
- MaxInRange(fMult, 0.2, mean, rms);
+ }
+ //__________________________________________________________________
+ /**
+ * Draw the multiplicity distribution
+ *
+ */
+ void DrawMult(TCanvas* /* c */, TLegend* leg)
+ {
+ fMult->GetXaxis()->SetRangeUser(0.2,20);
+ Double_t integral = fMult->Integral();
+ Info("DrawMult", "Integral in range [0.2,20]=%f (%f entries)",
+ integral, fMult->GetEntries());
+ fMult->Scale(1. / integral);
+ Info("DrawMult", "Integral in range [0.2,20]=%f (%f entries)",
+ fMult->Integral(), fMult->GetEntries());
+ Double_t max = 1.5 * fMult->GetMaximum();
+ fMult->GetXaxis()->SetRangeUser(0,4);
+ fMult->SetMaximum(max);
+ fMult->SetStats(kFALSE);
+ fMult->Draw("e same");
+ leg->AddEntry(fMult, "Strip signal", "lf");
+
+ }
+ //__________________________________________________________________
+ /**
+ * Find the minimum and maximum values in range
+ *
+ * @param xmin
+ * @param xmax
+ * @param ymax
+ */
+ void FindMinMax(Double_t& xmin, Double_t& xmax, Double_t& ymax)
+ {
+ fMult->GetXaxis()->SetRangeUser(0.1,4);
+ TSpectrum s(10);
+ Int_t nPeak = s.Search(fMult);
+ fMult->GetXaxis()->SetRangeUser(0.1, 4);
+ Int_t bmax = fMult->GetMaximumBin();
+ xmax = fMult->GetBinCenter(bmax);
+ fMult->GetXaxis()->SetRangeUser(0.1, xmax);
+ Double_t bmin = fMult->GetMinimumBin();
+ xmin = fMult->GetBinCenter(bmin);
+ ymax = fMult->GetBinContent(bmax);
+ Info("Finish", "Simple peak finding found x_max=%f, x_min=%f y_max=%g",
+ xmax, xmin, ymax);
+
+ if (nPeak > 1) {
+ TPolyMarker* pm = static_cast<TPolyMarker*>(fMult->GetListOfFunctions()
+ ->FindObject("TPolyMarker"));
+
+ // Peaks are ordered by size
+ Double_t* peakX = pm->GetX();
+ Double_t* peakY = pm->GetY();
+ Double_t xlmax = peakX[1];
+ xmax = peakX[0];
+ ymax = peakY[0];
+ if (xmax < xlmax) {
+ xmax = xlmax;
+ xlmax = peakX[0];
+ ymax = peakY[1];
+ }
+
+ fMult->GetXaxis()->SetRangeUser(xlmax, xmax);
+ bmin = fMult->GetMinimumBin();
+ xmin = fMult->GetBinCenter(bmin);
+ Info("Finish", "Spectrum peak finding found x_max=%f, x_min=%f y_max=%g",
+ xmax, xmin, ymax);
+ }
+ }
+ //__________________________________________________________________
+ /**
+ * Fit a landau and a landau+gaussian to the data
+ *
+ * @param xmin Minimum of peak range
+ * @param xmax Maximum of peak range
+ * @param ymax Y value in MIP peak
+ * @param leg Legend
+ */
+ void FitLandau(Double_t xmin, Double_t xmax, Double_t& ymax, TLegend* leg)
+ {
+ fMult->GetXaxis()->SetRangeUser(xmin, xmax+(xmax-xmin));
+ Double_t mean = fMult->GetMean();
+ Double_t var = fMult->GetRMS();
+ Info("Finish", "Peak range [%f,%f] mean=%f, var=%f",
+ xmin, xmax+(xmax-xmin), mean, var);
+ Double_t lowCut = mean-var;
+ Double_t hiCut = 4; // 2*mean;
+ Info("Finish", "Low cut set to %f", lowCut);
+ fMult->GetXaxis()->SetRangeUser(0, hiCut);
+
+ TF1* pl = MakeLandau(lowCut, hiCut, xmax, var/2);
+ fMult->Fit(pl, "NEM", "", lowCut, hiCut);
+ ymax = pl->GetMaximum();
+ Info("Finish", "Maximum of landau is at %f (y_max=%g)",
+ pl->GetMaximumX(), ymax);
+
+ TF1* gl = MakeFoldLandau(lowCut, hiCut, pl, xmax, var/2);
+ gl->SetLineColor(kRed+1);
+ // gl->SetParLimits(1, xmax-var, xmax+var);
+ fMult->Fit(gl, "NEM", "", lowCut, hiCut);
+ TF1* l = MakeLandau(lowCut, hiCut, xmax, var/2);
+ l->SetLineColor(kGreen+1);
+ l->SetParameters(gl->GetParameter(0),
+ gl->GetParameter(1),
+ gl->GetParameter(2));
+ TF1* g = new TF1("g", "gaus", lowCut, hiCut);
+ g->SetParNames("A", "#mu", "#sigma");
+ g->SetLineColor(kBlue+1);
+ g->SetParameters(gl->GetParameter(3)*gl->GetParameter(0),
+ gl->GetParameter(1),
+ gl->GetParameter(2));
+ fMult->GetXaxis()->SetRangeUser(0,4);
+ fMult->DrawCopy("E");
+ fMult->DrawCopy("HIST SAME");
+ pl->Draw("same");
+ gl->Draw("same");
+ g->Draw("Same");
+ l->Draw("Same");
+ fCleanup.Add(pl);
+ fCleanup.Add(l);
+ fCleanup.Add(g);
+ fCleanup.Add(gl);
+
+ ShowFit(.5, .90, "Landau", pl);
+ ShowFit(.5, .65, "Landau+Gaussian", gl);
+
+ leg->AddEntry(pl, Form("Landau fit - #chi^{2}/NDF=%f",
+ pl->GetChisquare()/pl->GetNDF()), "l");
+ leg->AddEntry(gl, Form("Landau+Gaussian fit - #chi^{2}/NDF=%f",
+ gl->GetChisquare()/gl->GetNDF()), "l");
+ leg->AddEntry(l, "Landau part", "l");
+ leg->AddEntry(g, "Gaussian part", "l");
+ }
+
+ //__________________________________________________________________
+ /**
+ * Superimpose the response graph from the RPP
+ *
+ * @param ymax Y value of multiplicity spectra in the landau peak
+ * @param leg Legend
+ */
+ void DrawResponse(Double_t xmax, Double_t ymax, TLegend* leg)
+ {
+ TGraph* resp = GetResp();
+ // TGraph* corr = GetCorr();
+ Double_t* x = resp->GetX();
+ Double_t* y = resp->GetY(); // [MeV cm^2/g]
+ TGraph* gr = new TGraph(resp->GetN());
+ gr->SetName(Form("%sCorr", resp->GetName()));
+ gr->SetTitle(resp->GetTitle());
+ gr->SetLineStyle(resp->GetLineStyle());
+ gr->SetLineColor(kMagenta+1);
+ gr->SetLineWidth(2);
+ TGraph* gr2 = new TGraph(resp->GetN());
+ gr2->SetName(Form("%sCorr", resp->GetName()));
+ gr2->SetTitle(resp->GetTitle());
+ gr2->SetLineStyle(resp->GetLineStyle());
+ gr2->SetLineColor(kCyan+1);
+ gr2->SetLineWidth(2);
+ // const Double_t rho = 2.33; // [g/cm^3]
+ // const Double_t thk = 0.320; // [cm]
+ const Double_t mip = 1.664; // [MeV cm^2/g]
+ // const Double_t bgm = 3.4601; // beta*gamma of a MIP
+ // Double_t xs2 = corr->Eval(bgm); // [1]
+ // Double_t xss = 1.1;
+ Double_t xs = 1/mip;
+ Double_t xxmax = 0;
+ Double_t yymax = 0;
+ for (Int_t i = 0; i < gr->GetN(); i++) {
+ if (y[i] > yymax) {
+ yymax = y[i];
+ xxmax = xs * x[i];
+ }
+ gr->SetPoint(i, x[i] * xs, y[i] * ymax);
+ }
+ Info("DrawResponse", "Maximum at x=%f (xmax=%f)", xxmax, xmax);
+ Double_t xs2 = xmax / xxmax;
+ Info("DrawResponse", "Correction factor: %f", xs2);
+ for (Int_t i = 0; i < gr->GetN(); i++) {
+ gr2->SetPoint(i, x[i] * xs * xs2, y[i] * ymax);
+ }
+ gr->Draw("C same");
+ gr2->Draw("C same");
+
+ leg->AddEntry(gr, "Response", "l");
+ leg->AddEntry(gr2, "Response", "l");
+ }
+ //__________________________________________________________________
+ /**
+ * Fit sum of landaus to the multiplicity distribution
+ *
+ * @param xmin Minimum of range
+ * @param xmax Maximum of range
+ * @param leg Legend
+ *
+ * @return Fitted function
+ */
+ TF1* FitMultiLandau(Double_t xmin, Double_t xmax, TLegend* leg)
+ {
+ fMult->GetXaxis()->SetRangeUser(xmin, xmax+(xmax-xmin));
+ Double_t mean = fMult->GetMean();
+ Double_t rms = fMult->GetRMS();
+
Double_t x1 = mean-rms/2; // .75; // .8; // .65 / fCorr;
Double_t x2 = mean+rms/2; // 1.3; // 1.7; // fCorr;
TF1* l1 = FitPeak(1, fMult, x1, x2);
std::cout << "Range: " << x1 << "-" << x3 << std::endl;
fMult->GetXaxis()->SetRangeUser(0, 4);
- fMult->Fit(f, "0QR", "", x1, x3);
- fMult->Fit(f, "ME0R", "E1", x1, x3);
- fMult->DrawClone("same hist");
+ fMult->Fit(f, "NQR", "", x1, x3);
+ fMult->Fit(f, "NMER", "E1", x1, x3);
l1->SetLineColor(3);
l1->SetLineWidth(2);
f->SetRange(0, 4);
f->Draw("same");
- TLegend* l = new TLegend(0.6, 0.6, .89, .89);
- l->AddEntry(l1, "1 particle Landau", "l");
- if (l2) l->AddEntry(l2, "2 particle Landau", "l");
- l->AddEntry(f, "1+2 particle Landau", "l");
- l->SetFillColor(0);
- l->Draw("same");
+ fCleanup.Add(l1);
+ if (l2) fCleanup.Add(l2);
+ fCleanup.Add(f);
-#if 0
- c = new TCanvas("c2", "Landaus");
+ leg->AddEntry(l1, "1 particle Landau", "l");
+ if (l2) leg->AddEntry(l2, "2 particle Landau", "l");
+ leg->AddEntry(f, "1+2 particle Landau", "l");
+
+ return f;
+ }
+ //__________________________________________________________________
+ /**
+ * Draw landau functions in a separate canvas
+ *
+ * @param f Multi-landau function
+ */
+ void DrawLandaus(TF1* f)
+ {
+ TCanvas* c = new TCanvas("c2", "Landaus");
// c->SetLogy();
fMult->DrawClone("axis");
f->Draw("same");
ll2->SetLineColor(4);
ll2->Draw("same");
- Double_t y1 = fMult->GetMinimum() * 1.1;
- Double_t y2 = fMult->GetMaximum() * .9;
- Double_t xc1 = p1[1]-3*p1[2];
- Double_t xc2 = p2[1]-2*p2[2];
- Double_t xc3 = p2[1]-2*p2[2]+diff;
+ Double_t diff = ll2->GetParameter(1)-ll1->GetParameter(1);
+ Double_t y1 = fMult->GetMinimum() * 1.1;
+ Double_t y2 = fMult->GetMaximum() * .9;
+ Double_t xc1 = p1[1]-3*p1[2];
+ Double_t xc2 = p2[1]-2*p2[2];
+ Double_t xc3 = p2[1]-2*p2[2]+diff;
TLine* c1 = new TLine(xc1, y1, xc1, y2);
c1->Draw("same");
TLine* c2 = new TLine(xc2, y1, xc2, y2);
TLine* c3 = new TLine(xc3, y1, xc3, y2);
c3->Draw("same");
- l = new TLegend(0.6, 0.6, .89, .89);
+ TLegend* l = new TLegend(0.6, 0.6, .89, .89);
l->AddEntry(ll1, "1 particle Landau", "l");
l->AddEntry(ll2, "2 particle Landau", "l");
l->AddEntry(f, "1+2 particle Landau", "l");
l->SetFillColor(0);
l->Draw("same");
-#endif
+ c->Modified();
+ c->Update();
+ c->cd();
+ }
- return kTRUE;
+ //__________________________________________________________________
+ /**
+ * Get the response functin @f$ f(\Delta_p/x)@f$ from Review of
+ * Particle Physics (fig. 27.7). It is scaled to the value at MPV.
+ *
+ * @return Graph of response
+ */
+ TGraph* GetResp()
+ {
+ static TGraph* graph = 0;
+ if (!graph) {
+ graph = new TGraph;
+ graph->SetName("si_resp");
+ graph->SetTitle("f(#Delta/x) scaled to the MPV value ");
+ graph->GetHistogram()->SetXTitle("#Delta/x (MeVcm^{2}/g)");
+ graph->GetHistogram()->SetYTitle("f(#Delta/x)");
+ graph->SetLineColor(kBlue+1);
+ graph->SetLineWidth(2);
+ graph->SetFillColor(kBlue+1);
+ graph->SetMarkerStyle(21);
+ graph->SetMarkerSize(0.6);
+#if 0
+ // Figure 27.7 or Review of Particle physics - Straggeling function in
+ // silicon of 500 MeV Pions, normalised to unity at the most probable
+ // value.
+ graph->SetPoint(0,0.808094,0.00377358);
+ graph->SetPoint(1,0.860313,0.0566038);
+ graph->SetPoint(2,0.891645,0.116981);
+ graph->SetPoint(3,0.912533,0.181132);
+ graph->SetPoint(4,0.928198,0.260377);
+ graph->SetPoint(5,0.938642,0.320755);
+ graph->SetPoint(6,0.954308,0.377358);
+ graph->SetPoint(7,0.964752,0.433962);
+ graph->SetPoint(8,0.975196,0.490566);
+ graph->SetPoint(9,0.98564,0.550943);
+ graph->SetPoint(10,0.996084,0.611321);
+ graph->SetPoint(11,1.00653,0.667925);
+ graph->SetPoint(12,1.02219,0.732075);
+ graph->SetPoint(13,1.03264,0.784906);
+ graph->SetPoint(14,1.0483,0.845283);
+ graph->SetPoint(15,1.06397,0.901887);
+ graph->SetPoint(16,1.09008,0.958491);
+ graph->SetPoint(17,1.10574,0.984906);
+ graph->SetPoint(18,1.13708,1);
+ graph->SetPoint(19,1.13708,1);
+ graph->SetPoint(20,1.15796,0.988679);
+ graph->SetPoint(21,1.17363,0.966038);
+ graph->SetPoint(22,1.19974,0.916981);
+ graph->SetPoint(23,1.2154,0.89434);
+ graph->SetPoint(24,1.23629,0.837736);
+ graph->SetPoint(25,1.2624,0.784906);
+ graph->SetPoint(26,1.28329,0.724528);
+ graph->SetPoint(27,1.3094,0.664151);
+ graph->SetPoint(28,1.32507,0.611321);
+ graph->SetPoint(29,1.3564,0.550943);
+ graph->SetPoint(30,1.41384,0.445283);
+ graph->SetPoint(31,1.44517,0.392453);
+ graph->SetPoint(32,1.48695,0.335849);
+ graph->SetPoint(33,1.52872,0.286792);
+ graph->SetPoint(34,1.58094,0.237736);
+ graph->SetPoint(35,1.63838,0.196226);
+ graph->SetPoint(36,1.68016,0.169811);
+ graph->SetPoint(37,1.75326,0.135849);
+ graph->SetPoint(38,1.81593,0.113208);
+ graph->SetPoint(39,1.89426,0.0981132);
+ graph->SetPoint(40,1.96214,0.0830189);
+ graph->SetPoint(41,2.0718,0.0641509);
+ graph->SetPoint(42,2.19191,0.0490566);
+ graph->SetPoint(43,2.31723,0.0415094);
+ graph->SetPoint(44,2.453,0.0301887);
+ graph->SetPoint(45,2.53133,0.0264151);
+ graph->SetPoint(46,2.57833,0.0264151);
+#else
+ graph->SetPoint(0,0.8115124,0.009771987);
+ graph->SetPoint(1,0.9198646,0.228013);
+ graph->SetPoint(2,0.996614,0.5895765);
+ graph->SetPoint(3,1.041761,0.8241042);
+ graph->SetPoint(4,1.059819,0.8794788);
+ graph->SetPoint(5,1.077878,0.9348534);
+ graph->SetPoint(6,1.100451,0.980456);
+ graph->SetPoint(7,1.141084,0.9967427);
+ graph->SetPoint(8,1.204289,0.9153094);
+ graph->SetPoint(9,1.276524,0.742671);
+ graph->SetPoint(10,1.402935,0.465798);
+ graph->SetPoint(11,1.515801,0.3029316);
+ graph->SetPoint(12,1.73702,0.1465798);
+ graph->SetPoint(13,1.985327,0.08143322);
+ graph->SetPoint(14,2.301354,0.04234528);
+ graph->SetPoint(15,2.56772,0.02931596);
+#endif
+ }
+ return graph;
+ }
+ //__________________________________________________________________
+ /**
+ * Get the correction to Bethe-Bloc from Review of Particle Physics
+ * (fig 27.8).
+ *
+ * @return correction graph
+ */
+ TGraph* GetCorr()
+ {
+ static TGraph* graph = 0;
+ if (!graph) {
+ graph = new TGraph(14);
+ graph->SetName("graph");
+ graph->SetTitle("(#Delta_{p}/x)/(dE/dx)|_{mip} for 320#mu Si");
+ graph->GetHistogram()->SetXTitle("#beta#gamma = p/m");
+ graph->GetHistogram()->SetYTitle("#frac{#Delta_{p}/x)}{(dE/dx)|_{mip}}");
+ graph->SetFillColor(1);
+ graph->SetLineColor(7);
+ graph->SetMarkerStyle(21);
+ graph->SetMarkerSize(0.6);
+ graph->SetPoint(0,1.196058,0.9944915);
+ graph->SetPoint(1,1.28502,0.9411017);
+ graph->SetPoint(2,1.484334,0.8559322);
+ graph->SetPoint(3,1.984617,0.7491525);
+ graph->SetPoint(4,2.658367,0.6983051);
+ graph->SetPoint(5,3.780227,0.6779661);
+ graph->SetPoint(6,4.997358,0.6741525);
+ graph->SetPoint(7,8.611026,0.684322);
+ graph->SetPoint(8,15.28296,0.6995763);
+ graph->SetPoint(9,41.54516,0.7186441);
+ graph->SetPoint(10,98.91461,0.7288136);
+ graph->SetPoint(11,203.2734,0.7326271);
+ graph->SetPoint(12,505.6421,0.7338983);
+ graph->SetPoint(13,896.973,0.7338983);
+ }
+ return graph;
+ }
+ //__________________________________________________________________
+ /**
+ * Make a Landau function object.
+ *
+ * @param min Minimum of fit range
+ * @param max Maximum of fit range
+ * @param p Peak position
+ * @param v Variance around peak
+ *
+ * @return Landau function object
+ */
+ TF1* MakeLandau(Double_t min, Double_t max, Double_t p, Double_t v)
+ {
+ TF1* f = new TF1("l", "landau", min, max);
+ f->SetParNames("A", "#delta", "#xi");
+ f->SetParameters(1, p, p/10);
+ if (false) f->FixParameter(1,p);
+ else if (false) f->SetParLimits(1, p-v, p+v);
+ return f;
+ }
+ //__________________________________________________________________
+ /**
+ * Make a Landau, folded with a gaussian, function object
+ *
+ * @param min Minimum of fit range
+ * @param max Maximum of fit range
+ * @param l Seed Landau function object
+ * @param p Peak position
+ * @param v Variance around peak
+ *
+ * @return Landau+Gaus function object
+ */
+ TF1* MakeFoldLandau(Double_t min, Double_t max, TF1* l,
+ Double_t p, Double_t v)
+ {
+ TF1* f = new TF1("gl", &foldLandau, min, max, 5);
+ f->SetParNames("A", "#delta", "#xi", "B", "#sigma");
+ f->SetParameters(l->GetParameter(0),
+ l->GetParameter(1),
+ l->GetParameter(2),
+ l->GetParameter(0)/1000,
+ l->GetParameter(2));
+ f->SetParLimits(3, 1e-7, 1e-1);
+ if (false) f->FixParameter(1,p);
+ else if (true) f->SetParLimits(1, p-2*v, p+2*v);
+ return f;
}
ClassDef(DrawESD,0);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff