Script to get # of dead channels from OCDB
[u/mrichter/AliRoot.git] / FMD / scripts / DrawESD.C
CommitLineData
2b893216 1//____________________________________________________________________
2//
3// $Id$
4//
5// Script that contains a class to draw eloss from hits, versus ADC
6// counts from digits, using the AliFMDInputHits class in the util library.
7//
8// It draws the energy loss versus the p/(mq^2). It can be overlayed
9// with the Bethe-Bloc curve to show how the simulation behaves
10// relative to the expected.
11//
12// Use the script `Compile.C' to compile this class using ACLic.
13//
14#include <TH1D.h>
15#include <AliFMDHit.h>
16#include <AliFMDDigit.h>
17#include <AliFMDInput.h>
18#include <AliFMDUShortMap.h>
19#include <AliFMDFloatMap.h>
20#include <AliFMDRecPoint.h>
21#include <AliESDFMD.h>
22#include <AliLog.h>
23#include <iostream>
24#include <TStyle.h>
25#include <TArrayF.h>
26#include <TCanvas.h>
27#include <TMath.h>
28#include <TF1.h>
29#include <TSpectrum.h>
30#include <TLegend.h>
2180cab3 31#include <TLatex.h>
2b893216 32#include <TLine.h>
2180cab3 33#include <TPolyMarker.h>
34#include <TSpectrum.h>
35#include <TList.h>
36#include <TGraph.h>
37
38Double_t landau(Double_t* xp, Double_t* pp)
39{
40 Double_t x = xp[0];
41 Double_t A = pp[0];
42 Double_t mpv = pp[1];
43 Double_t w = pp[2];
44
45 Double_t v = mpv; // + w * 0.22278;
46 return A * TMath::Landau(x, v, w);
47}
48
49Double_t foldLandau(Double_t* xp, Double_t* pp)
50{
51 Double_t x = xp[0];
52 Double_t A = pp[0];
53 Double_t mpv = pp[1];
54 Double_t w = pp[2];
55 Double_t B = pp[3];
56 Double_t sigma = pp[4];
57
58 return A * (TMath::Landau(x,mpv,w) + B * TMath::Gaus(x, mpv, sigma));
59}
2b893216 60
61/** @class DrawESD
62 @brief Draw digit ADC versus Rec point mult
63 @code
64 Root> .L Compile.C
65 Root> Compile("DrawESD.C")
66 Root> DrawESD c
67 Root> c.Run();
68 @endcode
69 @ingroup FMD_script
70 */
71class DrawESD : public AliFMDInput
72{
73private:
74 TH1D* fMult; // Histogram
75 const Double_t fCorr;
2180cab3 76 TList fCleanup;
2b893216 77public:
78 //__________________________________________________________________
2180cab3 79 DrawESD(Int_t n=2000, Double_t mmin=-0.5, Double_t mmax=15.5)
2b893216 80 : fCorr(1) // 0.68377 / 1.1)
81 {
82 AddLoad(kESD);
2180cab3 83 fMult = new TH1D("mult", "#DeltaE/#DeltaE_{MIP)", n, mmin, mmax);
69893a66 84 fMult->Sumw2();
2180cab3 85 fMult->SetFillColor(kRed+1);
86 fMult->SetFillStyle(3001);
87 fMult->SetXTitle("#DeltaE/#DeltaE_{MIP}");
88 fCleanup.Add(fMult);
89 }
90 ~DrawESD()
91 {
92 fCleanup.Delete();
2b893216 93 }
94 //__________________________________________________________________
95 /** Begining of event
96 @param ev Event number
97 @return @c false on error */
98 Bool_t Begin(Int_t ev)
99 {
100 return AliFMDInput::Begin(ev);
101 }
102 //__________________________________________________________________
2180cab3 103 Bool_t ProcessESD(UShort_t /* det */,
104 Char_t /* rng */,
105 UShort_t /* sec */,
106 UShort_t /* str */,
107 Float_t eta,
2b893216 108 Float_t mult)
109 {
110 // Cache the energy loss
2180cab3 111 // if (mult > 0)
112 // Info("ProcessESD", "FMD%d%c[%2d,%3d]=%f", det, rng, sec, str, mult);
113 if (mult <= 0 || mult == AliESDFMD::kInvalidMult) return kTRUE;
114 Double_t x = mult;
115 if (!fESD->IsAngleCorrected()) {
116 Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
117 Double_t corr = TMath::Abs(TMath::Cos(theta));
118 Double_t cmult = corr * mult;
119 x = cmult;
120 }
121 if (x > 0.001) fMult->Fill(x);
122 return kTRUE;
123 }
124 //__________________________________________________________________
125 Bool_t ProcessESDs()
126 {
127 if (!AliFMDInput::ProcessESDs()) return kFALSE;
128 // Info("ProcessESDs", "ESD is %sangle corrected",
129 // fESD->IsAngleCorrected() ? "" : "not ");
2b893216 130 return kTRUE;
131 }
132 //__________________________________________________________________
f48d9b11 133 TF1* FitPeak(Int_t n, TH1D* hist, Double_t min, Double_t& max)
134 {
faf80567 135 if (TMath::Abs(max-min) < .25) return 0;
2180cab3 136 std::cout << "Fit peak in range " << min << " to " << max << std::endl;
f48d9b11 137 TF1* l = new TF1(Form("l%02d", n), "landau", min, max);
138 hist->GetXaxis()->SetRangeUser(0, 4);
139 hist->Fit(l, "0Q", "", min, max);
140 Double_t mpv = l->GetParameter(1);
141 Double_t empv = l->GetParError(1);
142 Double_t sigma = l->GetParameter(2);
143 l->SetRange(mpv-empv, mpv+3*sigma);
144 hist->Fit(l, "EMQ0", "", mpv-3*empv, mpv+3*sigma);
145 std::cout << "Peak # " << n << " [" << min << "," << max << "]\n"
146 << " MPV: " << l->GetParameter(1)
147 << " +/- " << l->GetParError(1)
148 << " Var: " << l->GetParameter(2)
149 << " +/- " << l->GetParError(2)
150 << " Chi^2/NDF: " << l->GetChisquare() / l->GetNDF()
151 << std::endl;
152 mpv = l->GetParameter(1);
153 sigma = l->GetParameter(2);
154 min = mpv - sigma * 2; // (n==1 ? 3 : 2);
155 max = mpv + sigma * 3;
156 // l->SetRange(min, max);
157 l->Draw("same");
158 return l;
159 }
160 //__________________________________________________________________
161 void MaxInRange(TH1D* hist, Double_t min, Double_t& mean, Double_t& var)
162 {
163 hist->GetXaxis()->SetRangeUser(min, 4);
164 mean = hist->GetMean();
165 var = hist->GetRMS();
166 }
167
168 //__________________________________________________________________
2180cab3 169 const char* PrettyFloat(float x)
170 {
171 if (x == 0) return Form("%9.4f", x);
172 float e = TMath::Floor(TMath::Log10(TMath::Abs(x)));
173 if (TMath::Abs(e) < 4) {
174 return Form("%9.4f", x);
175 }
176 float f = x / TMath::Power(10,e);
177 return Form("%4.2f#times10^{%d}", f, int(e));
178 }
179 //__________________________________________________________________
180 void ShowFit(Double_t x1, Double_t y1, const char* title,
181 TF1* f, Double_t dx=0, Double_t dy=0.05)
182 {
183 Double_t x = x1, y = y1;
184 TLatex* latex = new TLatex(x, y, title);
185 latex->SetTextFont(132);
186 latex->SetTextSize(0.8*dy);
187 latex->SetNDC();
188 latex->Draw();
189 x -= dx;
190 y -= dy;
191 const Double_t eqDx=0.1;
192 Double_t chi2 = f->GetChisquare();
193 Int_t ndf = f->GetNDF();
194 Double_t prob = f->GetProb();
195 latex->DrawLatex(x, y, "#chi^{2}/NDF");
196 latex->DrawLatex(x+eqDx, y, Form("= %7.4f/%3d=%5.2f (%7.3f%%)",
197 chi2, ndf, chi2/ndf, 100*prob));
198 Int_t n = f->GetNpar();
199 Double_t* p = f->GetParameters();
200 Double_t* e = f->GetParErrors();
201 for (int i = 0; i < n; i++) {
202 x -= dx;
203 y -= dy;
204 latex->DrawLatex(x, y, f->GetParName(i));
205 latex->DrawLatex(x+eqDx, y, Form("= %s", PrettyFloat(p[i])));
206 latex->DrawLatex(x+2.2*eqDx, y, Form("#pm %s", PrettyFloat(e[i])));
207 }
208 }
209 //__________________________________________________________________
2b893216 210 Bool_t Finish()
211 {
faf80567 212 Info("Finish", "Will draw results");
2b893216 213 gStyle->SetPalette(1);
214 gStyle->SetOptTitle(0);
215 gStyle->SetOptFit(1111111);
216 gStyle->SetCanvasColor(0);
217 gStyle->SetCanvasBorderSize(0);
218 gStyle->SetPadColor(0);
219 gStyle->SetPadBorderSize(0);
220
faf80567 221 if (fMult->GetEntries() <= 0) return kFALSE;
222
2180cab3 223 TCanvas* c = new TCanvas("c", "C", 1200, 800);
224 fCleanup.Add(c);
2b893216 225 c->cd();
2180cab3 226 c->SetLogy();
227 c->SetTopMargin(0.05);
228 c->SetRightMargin(0.05);
229 c->SetFillColor(0);
230 c->SetBorderMode(0);
231
232 TLegend* leg = new TLegend(.1, .1, .4, .2);
233 leg->SetFillColor(0);
234 leg->SetBorderSize(1);
235
236 DrawMult(c, leg);
237
238 Double_t xmax=0, xmin=0, ymax=0;
239 FindMinMax(xmin, xmax, ymax);
240
241 FitLandau(xmin, xmax, ymax, leg);
242 DrawResponse(xmax, ymax, leg);
faf80567 243
2180cab3 244 // TF1* f = FitMultiLandau(xmin, xmax, leg);
245 // DrawLandaus(f);
246
247 c->cd();
248 leg->Draw();
249
250 c->Modified();
251 c->Update();
252 c->cd();
253 c->SaveAs("esd_eloss.png");
254
255 fCleanup.Add(leg);
256
faf80567 257 return kTRUE;
258
2180cab3 259 }
260 //__________________________________________________________________
261 /**
262 * Draw the multiplicity distribution
263 *
264 */
265 void DrawMult(TCanvas* /* c */, TLegend* leg)
266 {
267 fMult->GetXaxis()->SetRangeUser(0.2,20);
268 Double_t integral = fMult->Integral();
269 Info("DrawMult", "Integral in range [0.2,20]=%f (%f entries)",
270 integral, fMult->GetEntries());
271 fMult->Scale(1. / integral);
272 Info("DrawMult", "Integral in range [0.2,20]=%f (%f entries)",
273 fMult->Integral(), fMult->GetEntries());
274 Double_t max = 1.5 * fMult->GetMaximum();
275 fMult->GetXaxis()->SetRangeUser(0,4);
276 fMult->SetMaximum(max);
277 fMult->SetStats(kFALSE);
278 fMult->Draw("e same");
279 leg->AddEntry(fMult, "Strip signal", "lf");
280
281 }
282 //__________________________________________________________________
283 /**
284 * Find the minimum and maximum values in range
285 *
286 * @param xmin
287 * @param xmax
288 * @param ymax
289 */
290 void FindMinMax(Double_t& xmin, Double_t& xmax, Double_t& ymax)
291 {
292 fMult->GetXaxis()->SetRangeUser(0.1,4);
293 TSpectrum s(10);
294 Int_t nPeak = s.Search(fMult);
295 fMult->GetXaxis()->SetRangeUser(0.1, 4);
296 Int_t bmax = fMult->GetMaximumBin();
297 xmax = fMult->GetBinCenter(bmax);
298 fMult->GetXaxis()->SetRangeUser(0.1, xmax);
299 Double_t bmin = fMult->GetMinimumBin();
300 xmin = fMult->GetBinCenter(bmin);
301 ymax = fMult->GetBinContent(bmax);
302 Info("Finish", "Simple peak finding found x_max=%f, x_min=%f y_max=%g",
303 xmax, xmin, ymax);
304
305 if (nPeak > 1) {
306 TPolyMarker* pm = static_cast<TPolyMarker*>(fMult->GetListOfFunctions()
307 ->FindObject("TPolyMarker"));
308
309 // Peaks are ordered by size
310 Double_t* peakX = pm->GetX();
311 Double_t* peakY = pm->GetY();
312 Double_t xlmax = peakX[1];
313 xmax = peakX[0];
314 ymax = peakY[0];
315 if (xmax < xlmax) {
316 xmax = xlmax;
317 xlmax = peakX[0];
318 ymax = peakY[1];
319 }
320
321 fMult->GetXaxis()->SetRangeUser(xlmax, xmax);
322 bmin = fMult->GetMinimumBin();
323 xmin = fMult->GetBinCenter(bmin);
324 Info("Finish", "Spectrum peak finding found x_max=%f, x_min=%f y_max=%g",
325 xmax, xmin, ymax);
326 }
327 }
328 //__________________________________________________________________
329 /**
330 * Fit a landau and a landau+gaussian to the data
331 *
332 * @param xmin Minimum of peak range
333 * @param xmax Maximum of peak range
334 * @param ymax Y value in MIP peak
335 * @param leg Legend
336 */
337 void FitLandau(Double_t xmin, Double_t xmax, Double_t& ymax, TLegend* leg)
338 {
339 fMult->GetXaxis()->SetRangeUser(xmin, xmax+(xmax-xmin));
340 Double_t mean = fMult->GetMean();
341 Double_t var = fMult->GetRMS();
342 Info("Finish", "Peak range [%f,%f] mean=%f, var=%f",
343 xmin, xmax+(xmax-xmin), mean, var);
344 Double_t lowCut = mean-var;
345 Double_t hiCut = 4; // 2*mean;
346 Info("Finish", "Low cut set to %f", lowCut);
347 fMult->GetXaxis()->SetRangeUser(0, hiCut);
348
349 TF1* pl = MakeLandau(lowCut, hiCut, xmax, var/2);
350 fMult->Fit(pl, "NEM", "", lowCut, hiCut);
351 ymax = pl->GetMaximum();
352 Info("Finish", "Maximum of landau is at %f (y_max=%g)",
353 pl->GetMaximumX(), ymax);
354
355 TF1* gl = MakeFoldLandau(lowCut, hiCut, pl, xmax, var/2);
356 gl->SetLineColor(kRed+1);
357 // gl->SetParLimits(1, xmax-var, xmax+var);
358 fMult->Fit(gl, "NEM", "", lowCut, hiCut);
359 TF1* l = MakeLandau(lowCut, hiCut, xmax, var/2);
360 l->SetLineColor(kGreen+1);
361 l->SetParameters(gl->GetParameter(0),
362 gl->GetParameter(1),
363 gl->GetParameter(2));
364 TF1* g = new TF1("g", "gaus", lowCut, hiCut);
365 g->SetParNames("A", "#mu", "#sigma");
366 g->SetLineColor(kBlue+1);
367 g->SetParameters(gl->GetParameter(3)*gl->GetParameter(0),
368 gl->GetParameter(1),
369 gl->GetParameter(2));
370 fMult->GetXaxis()->SetRangeUser(0,4);
371 fMult->DrawCopy("E");
372 fMult->DrawCopy("HIST SAME");
373 pl->Draw("same");
374 gl->Draw("same");
375 g->Draw("Same");
376 l->Draw("Same");
377 fCleanup.Add(pl);
378 fCleanup.Add(l);
379 fCleanup.Add(g);
380 fCleanup.Add(gl);
381
382 ShowFit(.5, .90, "Landau", pl);
383 ShowFit(.5, .65, "Landau+Gaussian", gl);
384
385 leg->AddEntry(pl, Form("Landau fit - #chi^{2}/NDF=%f",
386 pl->GetChisquare()/pl->GetNDF()), "l");
387 leg->AddEntry(gl, Form("Landau+Gaussian fit - #chi^{2}/NDF=%f",
388 gl->GetChisquare()/gl->GetNDF()), "l");
389 leg->AddEntry(l, "Landau part", "l");
390 leg->AddEntry(g, "Gaussian part", "l");
391 }
392
393 //__________________________________________________________________
394 /**
395 * Superimpose the response graph from the RPP
396 *
397 * @param ymax Y value of multiplicity spectra in the landau peak
398 * @param leg Legend
399 */
400 void DrawResponse(Double_t xmax, Double_t ymax, TLegend* leg)
401 {
402 TGraph* resp = GetResp();
403 // TGraph* corr = GetCorr();
404 Double_t* x = resp->GetX();
405 Double_t* y = resp->GetY(); // [MeV cm^2/g]
406 TGraph* gr = new TGraph(resp->GetN());
407 gr->SetName(Form("%sCorr", resp->GetName()));
408 gr->SetTitle(resp->GetTitle());
409 gr->SetLineStyle(resp->GetLineStyle());
410 gr->SetLineColor(kMagenta+1);
411 gr->SetLineWidth(2);
412 TGraph* gr2 = new TGraph(resp->GetN());
413 gr2->SetName(Form("%sCorr", resp->GetName()));
414 gr2->SetTitle(resp->GetTitle());
415 gr2->SetLineStyle(resp->GetLineStyle());
416 gr2->SetLineColor(kCyan+1);
417 gr2->SetLineWidth(2);
418 // const Double_t rho = 2.33; // [g/cm^3]
419 // const Double_t thk = 0.320; // [cm]
420 const Double_t mip = 1.664; // [MeV cm^2/g]
421 // const Double_t bgm = 3.4601; // beta*gamma of a MIP
422 // Double_t xs2 = corr->Eval(bgm); // [1]
423 // Double_t xss = 1.1;
424 Double_t xs = 1/mip;
425 Double_t xxmax = 0;
426 Double_t yymax = 0;
427 for (Int_t i = 0; i < gr->GetN(); i++) {
428 if (y[i] > yymax) {
429 yymax = y[i];
430 xxmax = xs * x[i];
431 }
432 gr->SetPoint(i, x[i] * xs, y[i] * ymax);
433 }
434 Info("DrawResponse", "Maximum at x=%f (xmax=%f)", xxmax, xmax);
435 Double_t xs2 = xmax / xxmax;
436 Info("DrawResponse", "Correction factor: %f", xs2);
437 for (Int_t i = 0; i < gr->GetN(); i++) {
438 gr2->SetPoint(i, x[i] * xs * xs2, y[i] * ymax);
439 }
440 gr->Draw("C same");
441 gr2->Draw("C same");
442
443 leg->AddEntry(gr, "Response", "l");
444 leg->AddEntry(gr2, "Response", "l");
445 }
446 //__________________________________________________________________
447 /**
448 * Fit sum of landaus to the multiplicity distribution
449 *
450 * @param xmin Minimum of range
451 * @param xmax Maximum of range
452 * @param leg Legend
453 *
454 * @return Fitted function
455 */
456 TF1* FitMultiLandau(Double_t xmin, Double_t xmax, TLegend* leg)
457 {
458 fMult->GetXaxis()->SetRangeUser(xmin, xmax+(xmax-xmin));
459 Double_t mean = fMult->GetMean();
460 Double_t rms = fMult->GetRMS();
461
f48d9b11 462 Double_t x1 = mean-rms/2; // .75; // .8; // .65 / fCorr;
463 Double_t x2 = mean+rms/2; // 1.3; // 1.7; // fCorr;
464 TF1* l1 = FitPeak(1, fMult, x1, x2);
465 x2 = TMath::Max(mean+rms/2, x2);
466 MaxInRange(fMult, x2, mean, rms);
467 Double_t x3 = mean + rms;
468 TF1* l2 = FitPeak(2, fMult, x2, x3);
faf80567 469 TF1* f = 0;
470 Double_t diff = 0;
471 if (l2) {
472 diff = l2->GetParameter(1)-l1->GetParameter(1);
473 f = new TF1("user", "landau(0)+landau(3)", x1, x3);
474 f->SetParNames("A_{1}", "Mpv_{1}", "#sigma_{1}",
475 "A_{2}", "Mpv_{2}", "#sigma_{2}",
476 "A_{3}", "Mpv_{3}", "#sigma_{3}");
477 f->SetParameters(l1->GetParameter(0),
478 l1->GetParameter(1),
479 l1->GetParameter(2),
480 l2 ? l2->GetParameter(0) : 0,
481 l2 ? l2->GetParameter(1) : 0,
482 l2 ? l2->GetParameter(2) : 0,
483 l2->GetParameter(0)/10,
484 l2->GetParameter(1) + diff,
485 l2->GetParameter(2));
486 }
487 else {
488 x3 = x2;
489 f = new TF1("usr", "landau", x1, x3);
490 }
491
492 std::cout << "Range: " << x1 << "-" << x3 << std::endl;
493
2b893216 494 fMult->GetXaxis()->SetRangeUser(0, 4);
2180cab3 495 fMult->Fit(f, "NQR", "", x1, x3);
496 fMult->Fit(f, "NMER", "E1", x1, x3);
2b893216 497
2b893216 498 l1->SetLineColor(3);
499 l1->SetLineWidth(2);
500 l1->SetRange(0, 4);
501 l1->Draw("same");
faf80567 502 if (l2) {
503 l2->SetLineColor(4);
504 l2->SetLineWidth(2);
505 l2->SetRange(0, 4);
506 l2->Draw("same");
507 }
2b893216 508 f->SetLineWidth(2);
509 f->SetRange(0, 4);
510 f->Draw("same");
511
2180cab3 512 fCleanup.Add(l1);
513 if (l2) fCleanup.Add(l2);
514 fCleanup.Add(f);
2b893216 515
516
2180cab3 517 leg->AddEntry(l1, "1 particle Landau", "l");
518 if (l2) leg->AddEntry(l2, "2 particle Landau", "l");
519 leg->AddEntry(f, "1+2 particle Landau", "l");
520
521 return f;
522 }
523 //__________________________________________________________________
524 /**
525 * Draw landau functions in a separate canvas
526 *
527 * @param f Multi-landau function
528 */
529 void DrawLandaus(TF1* f)
530 {
531 TCanvas* c = new TCanvas("c2", "Landaus");
faf80567 532 // c->SetLogy();
2b893216 533 fMult->DrawClone("axis");
534 f->Draw("same");
f48d9b11 535 Double_t* p1 = f->GetParameters();
536 Double_t* p2 = &(p1[3]);
2b893216 537 TF1* ll1 = new TF1("ll1", "landau", 0, 4);
f48d9b11 538 ll1->SetParameters(p1);
2b893216 539 ll1->SetLineColor(3);
540 ll1->Draw("same");
541 TF1* ll2 = new TF1("ll2", "landau", 0, 4);
f48d9b11 542 ll2->SetParameters(p2);
2b893216 543 ll2->SetLineColor(4);
544 ll2->Draw("same");
545
2180cab3 546 Double_t diff = ll2->GetParameter(1)-ll1->GetParameter(1);
547 Double_t y1 = fMult->GetMinimum() * 1.1;
548 Double_t y2 = fMult->GetMaximum() * .9;
549 Double_t xc1 = p1[1]-3*p1[2];
550 Double_t xc2 = p2[1]-2*p2[2];
551 Double_t xc3 = p2[1]-2*p2[2]+diff;
2b893216 552 TLine* c1 = new TLine(xc1, y1, xc1, y2);
553 c1->Draw("same");
554 TLine* c2 = new TLine(xc2, y1, xc2, y2);
555 c2->Draw("same");
f48d9b11 556 TLine* c3 = new TLine(xc3, y1, xc3, y2);
557 c3->Draw("same");
2b893216 558
2180cab3 559 TLegend* l = new TLegend(0.6, 0.6, .89, .89);
2b893216 560 l->AddEntry(ll1, "1 particle Landau", "l");
561 l->AddEntry(ll2, "2 particle Landau", "l");
562 l->AddEntry(f, "1+2 particle Landau", "l");
563 l->SetFillColor(0);
564 l->Draw("same");
2b893216 565
2180cab3 566 c->Modified();
567 c->Update();
568 c->cd();
569 }
f48d9b11 570
2180cab3 571 //__________________________________________________________________
572 /**
573 * Get the response functin @f$ f(\Delta_p/x)@f$ from Review of
574 * Particle Physics (fig. 27.7). It is scaled to the value at MPV.
575 *
576 * @return Graph of response
577 */
578 TGraph* GetResp()
579 {
580 static TGraph* graph = 0;
581 if (!graph) {
582 graph = new TGraph;
583 graph->SetName("si_resp");
584 graph->SetTitle("f(#Delta/x) scaled to the MPV value ");
585 graph->GetHistogram()->SetXTitle("#Delta/x (MeVcm^{2}/g)");
586 graph->GetHistogram()->SetYTitle("f(#Delta/x)");
587 graph->SetLineColor(kBlue+1);
588 graph->SetLineWidth(2);
589 graph->SetFillColor(kBlue+1);
590 graph->SetMarkerStyle(21);
591 graph->SetMarkerSize(0.6);
592#if 0
593 // Figure 27.7 or Review of Particle physics - Straggeling function in
594 // silicon of 500 MeV Pions, normalised to unity at the most probable
595 // value.
596 graph->SetPoint(0,0.808094,0.00377358);
597 graph->SetPoint(1,0.860313,0.0566038);
598 graph->SetPoint(2,0.891645,0.116981);
599 graph->SetPoint(3,0.912533,0.181132);
600 graph->SetPoint(4,0.928198,0.260377);
601 graph->SetPoint(5,0.938642,0.320755);
602 graph->SetPoint(6,0.954308,0.377358);
603 graph->SetPoint(7,0.964752,0.433962);
604 graph->SetPoint(8,0.975196,0.490566);
605 graph->SetPoint(9,0.98564,0.550943);
606 graph->SetPoint(10,0.996084,0.611321);
607 graph->SetPoint(11,1.00653,0.667925);
608 graph->SetPoint(12,1.02219,0.732075);
609 graph->SetPoint(13,1.03264,0.784906);
610 graph->SetPoint(14,1.0483,0.845283);
611 graph->SetPoint(15,1.06397,0.901887);
612 graph->SetPoint(16,1.09008,0.958491);
613 graph->SetPoint(17,1.10574,0.984906);
614 graph->SetPoint(18,1.13708,1);
615 graph->SetPoint(19,1.13708,1);
616 graph->SetPoint(20,1.15796,0.988679);
617 graph->SetPoint(21,1.17363,0.966038);
618 graph->SetPoint(22,1.19974,0.916981);
619 graph->SetPoint(23,1.2154,0.89434);
620 graph->SetPoint(24,1.23629,0.837736);
621 graph->SetPoint(25,1.2624,0.784906);
622 graph->SetPoint(26,1.28329,0.724528);
623 graph->SetPoint(27,1.3094,0.664151);
624 graph->SetPoint(28,1.32507,0.611321);
625 graph->SetPoint(29,1.3564,0.550943);
626 graph->SetPoint(30,1.41384,0.445283);
627 graph->SetPoint(31,1.44517,0.392453);
628 graph->SetPoint(32,1.48695,0.335849);
629 graph->SetPoint(33,1.52872,0.286792);
630 graph->SetPoint(34,1.58094,0.237736);
631 graph->SetPoint(35,1.63838,0.196226);
632 graph->SetPoint(36,1.68016,0.169811);
633 graph->SetPoint(37,1.75326,0.135849);
634 graph->SetPoint(38,1.81593,0.113208);
635 graph->SetPoint(39,1.89426,0.0981132);
636 graph->SetPoint(40,1.96214,0.0830189);
637 graph->SetPoint(41,2.0718,0.0641509);
638 graph->SetPoint(42,2.19191,0.0490566);
639 graph->SetPoint(43,2.31723,0.0415094);
640 graph->SetPoint(44,2.453,0.0301887);
641 graph->SetPoint(45,2.53133,0.0264151);
642 graph->SetPoint(46,2.57833,0.0264151);
643#else
644 graph->SetPoint(0,0.8115124,0.009771987);
645 graph->SetPoint(1,0.9198646,0.228013);
646 graph->SetPoint(2,0.996614,0.5895765);
647 graph->SetPoint(3,1.041761,0.8241042);
648 graph->SetPoint(4,1.059819,0.8794788);
649 graph->SetPoint(5,1.077878,0.9348534);
650 graph->SetPoint(6,1.100451,0.980456);
651 graph->SetPoint(7,1.141084,0.9967427);
652 graph->SetPoint(8,1.204289,0.9153094);
653 graph->SetPoint(9,1.276524,0.742671);
654 graph->SetPoint(10,1.402935,0.465798);
655 graph->SetPoint(11,1.515801,0.3029316);
656 graph->SetPoint(12,1.73702,0.1465798);
657 graph->SetPoint(13,1.985327,0.08143322);
658 graph->SetPoint(14,2.301354,0.04234528);
659 graph->SetPoint(15,2.56772,0.02931596);
660#endif
661 }
662 return graph;
663 }
664 //__________________________________________________________________
665 /**
666 * Get the correction to Bethe-Bloc from Review of Particle Physics
667 * (fig 27.8).
668 *
669 * @return correction graph
670 */
671 TGraph* GetCorr()
672 {
673 static TGraph* graph = 0;
674 if (!graph) {
675 graph = new TGraph(14);
676 graph->SetName("graph");
677 graph->SetTitle("(#Delta_{p}/x)/(dE/dx)|_{mip} for 320#mu Si");
678 graph->GetHistogram()->SetXTitle("#beta#gamma = p/m");
679 graph->GetHistogram()->SetYTitle("#frac{#Delta_{p}/x)}{(dE/dx)|_{mip}}");
680 graph->SetFillColor(1);
681 graph->SetLineColor(7);
682 graph->SetMarkerStyle(21);
683 graph->SetMarkerSize(0.6);
684 graph->SetPoint(0,1.196058,0.9944915);
685 graph->SetPoint(1,1.28502,0.9411017);
686 graph->SetPoint(2,1.484334,0.8559322);
687 graph->SetPoint(3,1.984617,0.7491525);
688 graph->SetPoint(4,2.658367,0.6983051);
689 graph->SetPoint(5,3.780227,0.6779661);
690 graph->SetPoint(6,4.997358,0.6741525);
691 graph->SetPoint(7,8.611026,0.684322);
692 graph->SetPoint(8,15.28296,0.6995763);
693 graph->SetPoint(9,41.54516,0.7186441);
694 graph->SetPoint(10,98.91461,0.7288136);
695 graph->SetPoint(11,203.2734,0.7326271);
696 graph->SetPoint(12,505.6421,0.7338983);
697 graph->SetPoint(13,896.973,0.7338983);
698 }
699 return graph;
700 }
701 //__________________________________________________________________
702 /**
703 * Make a Landau function object.
704 *
705 * @param min Minimum of fit range
706 * @param max Maximum of fit range
707 * @param p Peak position
708 * @param v Variance around peak
709 *
710 * @return Landau function object
711 */
712 TF1* MakeLandau(Double_t min, Double_t max, Double_t p, Double_t v)
713 {
714 TF1* f = new TF1("l", "landau", min, max);
715 f->SetParNames("A", "#delta", "#xi");
716 f->SetParameters(1, p, p/10);
717 if (false) f->FixParameter(1,p);
718 else if (false) f->SetParLimits(1, p-v, p+v);
719 return f;
720 }
721 //__________________________________________________________________
722 /**
723 * Make a Landau, folded with a gaussian, function object
724 *
725 * @param min Minimum of fit range
726 * @param max Maximum of fit range
727 * @param l Seed Landau function object
728 * @param p Peak position
729 * @param v Variance around peak
730 *
731 * @return Landau+Gaus function object
732 */
733 TF1* MakeFoldLandau(Double_t min, Double_t max, TF1* l,
734 Double_t p, Double_t v)
735 {
736 TF1* f = new TF1("gl", &foldLandau, min, max, 5);
737 f->SetParNames("A", "#delta", "#xi", "B", "#sigma");
738 f->SetParameters(l->GetParameter(0),
739 l->GetParameter(1),
740 l->GetParameter(2),
741 l->GetParameter(0)/1000,
742 l->GetParameter(2));
743 f->SetParLimits(3, 1e-7, 1e-1);
744 if (false) f->FixParameter(1,p);
745 else if (true) f->SetParLimits(1, p-2*v, p+2*v);
746 return f;
2b893216 747 }
748
749 ClassDef(DrawESD,0);
750
751};
752
753//____________________________________________________________________
754//
755// EOF
756//