//____________________________________________________________________ // // $Id$ // // Script that contains a class to draw eloss from hits, versus ADC // counts from digits, using the AliFMDInputHits class in the util library. // // It draws the energy loss versus the p/(mq^2). It can be overlayed // with the Bethe-Bloc curve to show how the simulation behaves // relative to the expected. // // Use the script `Compile.C' to compile this class using ACLic. // #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include 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 @code Root> .L Compile.C Root> Compile("DrawESD.C") Root> DrawESD c Root> c.Run(); @endcode @ingroup FMD_script */ class DrawESD : public AliFMDInput { private: TH1D* fMult; // Histogram const Double_t fCorr; TList fCleanup; public: //__________________________________________________________________ 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", "#DeltaE/#DeltaE_{MIP)", n, mmin, mmax); fMult->Sumw2(); fMult->SetFillColor(kRed+1); fMult->SetFillStyle(3001); fMult->SetXTitle("#DeltaE/#DeltaE_{MIP}"); fCleanup.Add(fMult); } ~DrawESD() { fCleanup.Delete(); } //__________________________________________________________________ /** Begining of event @param ev Event number @return @c false on error */ Bool_t Begin(Int_t ev) { return AliFMDInput::Begin(ev); } //__________________________________________________________________ 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 <= 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 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); Double_t mpv = l->GetParameter(1); Double_t empv = l->GetParError(1); Double_t sigma = l->GetParameter(2); l->SetRange(mpv-empv, mpv+3*sigma); hist->Fit(l, "EMQ0", "", mpv-3*empv, mpv+3*sigma); std::cout << "Peak # " << n << " [" << min << "," << max << "]\n" << " MPV: " << l->GetParameter(1) << " +/- " << l->GetParError(1) << " Var: " << l->GetParameter(2) << " +/- " << l->GetParError(2) << " Chi^2/NDF: " << l->GetChisquare() / l->GetNDF() << std::endl; mpv = l->GetParameter(1); sigma = l->GetParameter(2); min = mpv - sigma * 2; // (n==1 ? 3 : 2); max = mpv + sigma * 3; // l->SetRange(min, max); l->Draw("same"); return l; } //__________________________________________________________________ void MaxInRange(TH1D* hist, Double_t min, Double_t& mean, Double_t& var) { hist->GetXaxis()->SetRangeUser(min, 4); mean = hist->GetMean(); 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() { Info("Finish", "Will draw results"); gStyle->SetPalette(1); gStyle->SetOptTitle(0); gStyle->SetOptFit(1111111); gStyle->SetCanvasColor(0); gStyle->SetCanvasBorderSize(0); gStyle->SetPadColor(0); gStyle->SetPadBorderSize(0); if (fMult->GetEntries() <= 0) return kFALSE; TCanvas* c = new TCanvas("c", "C", 1200, 800); fCleanup.Add(c); c->cd(); 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; } //__________________________________________________________________ /** * 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(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); x2 = TMath::Max(mean+rms/2, x2); MaxInRange(fMult, x2, mean, rms); Double_t x3 = mean + rms; TF1* l2 = FitPeak(2, fMult, x2, x3); TF1* f = 0; Double_t diff = 0; if (l2) { diff = l2->GetParameter(1)-l1->GetParameter(1); f = new TF1("user", "landau(0)+landau(3)", x1, x3); f->SetParNames("A_{1}", "Mpv_{1}", "#sigma_{1}", "A_{2}", "Mpv_{2}", "#sigma_{2}", "A_{3}", "Mpv_{3}", "#sigma_{3}"); f->SetParameters(l1->GetParameter(0), l1->GetParameter(1), l1->GetParameter(2), l2 ? l2->GetParameter(0) : 0, l2 ? l2->GetParameter(1) : 0, l2 ? l2->GetParameter(2) : 0, l2->GetParameter(0)/10, l2->GetParameter(1) + diff, l2->GetParameter(2)); } else { x3 = x2; f = new TF1("usr", "landau", x1, x3); } std::cout << "Range: " << x1 << "-" << x3 << std::endl; fMult->GetXaxis()->SetRangeUser(0, 4); fMult->Fit(f, "NQR", "", x1, x3); fMult->Fit(f, "NMER", "E1", x1, x3); l1->SetLineColor(3); l1->SetLineWidth(2); l1->SetRange(0, 4); l1->Draw("same"); if (l2) { l2->SetLineColor(4); l2->SetLineWidth(2); l2->SetRange(0, 4); l2->Draw("same"); } f->SetLineWidth(2); f->SetRange(0, 4); f->Draw("same"); fCleanup.Add(l1); if (l2) fCleanup.Add(l2); fCleanup.Add(f); 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"); Double_t* p1 = f->GetParameters(); Double_t* p2 = &(p1[3]); TF1* ll1 = new TF1("ll1", "landau", 0, 4); ll1->SetParameters(p1); ll1->SetLineColor(3); ll1->Draw("same"); TF1* ll2 = new TF1("ll2", "landau", 0, 4); ll2->SetParameters(p2); ll2->SetLineColor(4); ll2->Draw("same"); 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); c2->Draw("same"); TLine* c3 = new TLine(xc3, y1, xc3, y2); c3->Draw("same"); 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"); c->Modified(); c->Update(); c->cd(); } //__________________________________________________________________ /** * 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); }; //____________________________________________________________________ // // EOF //