Fixed some warning messages on request from FCA
[u/mrichter/AliRoot.git] / FMD / scripts / DrawESD.C
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>
31 #include <TLatex.h>
32 #include <TLine.h>
33 #include <TPolyMarker.h>
34 #include <TSpectrum.h>
35 #include <TList.h>
36 #include <TGraph.h>
37
38 Double_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
49 Double_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 }
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  */
71 class DrawESD : public AliFMDInput
72 {
73 private:
74   TH1D* fMult; // Histogram 
75   const Double_t fCorr;
76   TList fCleanup;
77 public:
78   //__________________________________________________________________
79   DrawESD(Int_t n=2000, Double_t mmin=-0.5, Double_t mmax=15.5) 
80     : fCorr(1) // 0.68377 / 1.1)
81   { 
82     AddLoad(kESD);
83     fMult = new TH1D("mult", "#DeltaE/#DeltaE_{MIP)", n, mmin, mmax);
84     fMult->Sumw2();
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();
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   //__________________________________________________________________
103   Bool_t ProcessESD(UShort_t /* det */, 
104                     Char_t   /* rng */, 
105                     UShort_t /* sec */, 
106                     UShort_t /* str */, 
107                     Float_t  eta, 
108                     Float_t  mult)
109   {
110     // Cache the energy loss 
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 ");
130     return kTRUE;
131   }
132   //__________________________________________________________________
133   TF1* FitPeak(Int_t n, TH1D* hist, Double_t min, Double_t& max)
134   {
135     if (TMath::Abs(max-min) < .25) return 0;
136     std::cout << "Fit peak in range " << min << " to " << max << std::endl;
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   //__________________________________________________________________
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   //__________________________________________________________________
210   Bool_t Finish()
211   {
212     Info("Finish", "Will draw results");
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     
221     if (fMult->GetEntries() <= 0) return kFALSE;
222     
223     TCanvas* c = new TCanvas("c", "C", 1200, 800);
224     fCleanup.Add(c);
225     c->cd();
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);
243     
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
257     return kTRUE;
258     
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
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);
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     
494     fMult->GetXaxis()->SetRangeUser(0, 4);
495     fMult->Fit(f, "NQR", "", x1, x3);
496     fMult->Fit(f, "NMER", "E1", x1, x3);
497
498     l1->SetLineColor(3);
499     l1->SetLineWidth(2);
500     l1->SetRange(0, 4);
501     l1->Draw("same");
502     if (l2) {
503       l2->SetLineColor(4);
504       l2->SetLineWidth(2);
505       l2->SetRange(0, 4);
506       l2->Draw("same");
507     }
508     f->SetLineWidth(2);
509     f->SetRange(0, 4);
510     f->Draw("same");
511
512     fCleanup.Add(l1);
513     if (l2) fCleanup.Add(l2);
514     fCleanup.Add(f);
515
516
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");
532     // c->SetLogy();
533     fMult->DrawClone("axis");
534     f->Draw("same");
535     Double_t* p1 = f->GetParameters();
536     Double_t* p2 = &(p1[3]);
537     TF1* ll1 = new TF1("ll1", "landau", 0, 4);
538     ll1->SetParameters(p1);
539     ll1->SetLineColor(3);
540     ll1->Draw("same");
541     TF1* ll2 = new TF1("ll2", "landau", 0, 4);
542     ll2->SetParameters(p2);
543     ll2->SetLineColor(4);
544     ll2->Draw("same");
545
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;
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");
556     TLine* c3 = new TLine(xc3, y1, xc3, y2);
557     c3->Draw("same");
558
559     TLegend* l = new TLegend(0.6, 0.6, .89, .89);
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");
565
566     c->Modified();
567     c->Update();
568     c->cd();
569   }
570
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;
747   }
748
749   ClassDef(DrawESD,0);
750   
751 };
752
753 //____________________________________________________________________
754 //
755 // EOF
756 //