10 static Double_t function(Double_t* x, Double_t* par)
12 const Double_t invSq2Pi = 1. / TMath::Sqrt(TMath::TwoPi());
13 const Double_t mpShift = -0.22278298;
15 const Double_t sc = 5;
20 Double_t mpc = par[1] - mpShift * xi;
21 Double_t area = par[2];
22 Double_t sigma = par[3];
23 Double_t xLow = x0 - sc * sigma;
24 Double_t xUpp = x0 + sc * sigma;
25 Double_t step = (xUpp - xLow) / np;
27 for (Int_t i = 1; i <= np/2; i++) {
28 Double_t x1 = xLow + (i - .5) * step;
29 Double_t x2 = xUpp - (i - .5) * step;
31 sum += TMath::Landau(x1, mpc, xi, true) * TMath::Gaus(x0,x1,sigma);
32 sum += TMath::Landau(x2, mpc, xi, true) * TMath::Gaus(x0,x2,sigma);
35 return area * step * sum * invSq2Pi / sigma;
38 static TF1* fit(TH1* hist,
48 TF1* func = new TF1(Form("func%s", hist->GetName()), &function,
49 range[0], range[1], 4);
50 func->SetParameters(guess);
51 func->SetParNames("#xi", "#Delta_{p}", "C", "#sigma");
52 for (Int_t i = 0; i < 4; i++)
53 func->SetParLimits(i, lowLimits[i], uppLimits[i]);
55 hist->Fit(func, "RB0");
57 for (Int_t i = 0; i < 4; i++) {
58 params[i] = func->GetParameter(i);
59 errors[i] = func->GetParError(i);
61 chi2 = func->GetChisquare();
66 static Bool_t search(TF1* f, Double_t start, Double_t step,
67 Bool_t peak, Double_t& res)
69 const Int_t maxCalls = 10000;
71 Double_t l = peak ? -1 : -1e300;
76 while ((i < maxCalls) &&
77 TMath::Abs(l - lOld) > 1e-6 &&
78 TMath::Abs(step) > 1e-8) {
83 if (!peak) l = TMath::Abs(l-res);
84 // Printf("mode=%d x=%g step=%g l=%g lOld=%g", mode, x, step, l, lOld);
86 if ((peak && l < lOld) || (!peak && l > lOld))
87 step = -step/10; // Go the other way in smaller steps
91 if (i >= maxCalls) return false;
95 static Int_t find(TF1* f, Int_t iXi, Int_t iMpc,
96 Double_t& maxX, Double_t& fwhm)
99 Double_t xi = f->GetParameter(iXi); // params[0];
100 Double_t mpc = f->GetParameter(iMpc); // params[1];
101 if (!search(f, mpc-0.1*xi, 0.05 * xi, true, maxX)) return false;
102 Double_t left = f->Eval(maxX) / 2;
103 Double_t right = left;
104 if (!search(f, maxX-.5*xi, xi, false, left)) return false;
105 if (!search(f, maxX+xi, -xi, false, right)) return false;
112 Int_t data[100] = {0, 0, 0, 0, 0, 0, 2, 6, 11, 18,
113 18, 55, 90,141,255,323,454,563,681,737,
114 821,796,832,720,637,558,519,460,357,291,
115 279,241,212,153,164,139,106, 95, 91, 76,
116 80, 80, 59, 58, 51, 30, 49, 23, 35, 28,
117 23, 22, 27, 27, 24, 20, 16, 17, 14, 20,
118 12, 12, 13, 10, 17, 7, 6, 12, 6, 12,
119 4, 9, 9, 10, 3, 4, 5, 2, 4, 1,
120 5, 5, 1, 7, 1, 6, 3, 3, 3, 4,
121 5, 4, 4, 2, 2, 7, 2, 4};
122 TH1D *ret = new TH1D("snr","Signal-to-noise",100,0,100);
123 for (Int_t i = 0; i < 100; i++) ret->Fill(i, data[i]);
125 ret->SetFillColor(kRed+1);
126 ret->SetFillStyle(3001);
127 ret->SetXTitle("#Delta");
128 ret->SetDirectory(0);
132 static void maxFwhm(TF1* f, Int_t iXi, Int_t iMpc)
134 printf("Find max and FWHM ...");
137 if (!find(f, iXi, iMpc, maxX, fwhm))
140 Printf(" Max: %f FWHM: %f", maxX, fwhm);
141 Double_t y = f->Eval(maxX);
142 TGraph* maxG = new TGraph(1);
143 maxG->SetPoint(0, maxX, y);
144 maxG->SetMarkerStyle(20);
145 maxG->SetMarkerColor(kBlue+2);
148 TLatex* maxT = new TLatex(1.2*maxX, y, Form("Max: %f", maxX));
149 maxT->SetTextAlign(11);
152 TGraph* fwhmG = new TGraph(2);
153 fwhmG->SetPoint(0, maxX-fwhm/2, y/2);
154 fwhmG->SetPoint(1, maxX+fwhm/2, y/2);
155 fwhmG->SetMarkerStyle(21);
156 fwhmG->SetMarkerColor(kMagenta+2);
159 TLatex* fwhmT = new TLatex(1.2*(maxX+fwhm/2), y/2,
160 Form("FWHM: %f", fwhm));
161 fwhmT->SetTextAlign(11);
170 Double_t range[] = { .3 * hist->GetMean(), 3 * hist->GetMean() };
171 Double_t guess[] = { 1.8, 20., 5e4, 3.0 };
172 Double_t low[] = { 0.5, 5., 1., 0.4 };
173 Double_t upp[] = { 5.0, 50., 1e6, 5.0 };
179 printf("Fitting ...");
180 TF1* f = fit(hist, range, guess, low, upp, params, errors, chi2, ndf);
184 printf("Drawing ...");
185 gStyle->SetOptStat(1111);
186 gStyle->SetOptFit(111);