6845cf4474f68450489cefdac2cbe1c7fe525c9b
[u/mrichter/AliRoot.git] / FMD / scripts / DrawHits.C
1 //____________________________________________________________________
2 //
3 // $Id$
4 //
5 // Script that contains a class to draw hits, using the
6 // 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 <TH2D.h>
15 #include <AliFMDHit.h>
16 #include <AliFMDInput.h>
17 #include <iostream>
18 #include <TStyle.h>
19 #include <TArrayF.h>
20 #include <TParticle.h>
21 #include <TCanvas.h>
22 #include <TGraphErrors.h>
23 #include <TDatabasePDG.h>
24 #include <TParticlePDG.h>
25 #include <TLegend.h>
26 #include <TArrow.h>
27 #include <TLatex.h>
28 #include <TF1.h>
29
30 /** @class DrawHits
31     @brief Draw hit energy loss
32     @code 
33     Root> .L Compile.C
34     Root> Compile("DrawHits.C")
35     Root> DrawHits c
36     Root> c.Run();
37     @endcode
38     @ingroup FMD_script
39  */
40 class DrawHits : public AliFMDInput
41 {
42 private:
43   TH2D* fElossVsPMQ; // Histogram 
44   TH1D* fEloss;
45   TParticlePDG* fPdg;
46   const Double_t fRho;
47   const Double_t fBetaGammaMip;
48 public:
49   //__________________________________________________________________
50   TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max) 
51   {
52     TArrayF bins(n+1);
53     bins[0]      = min;
54     if (n <= 20) {
55       for (Int_t i = 1; i < n+1; i++) bins[i] = bins[i-1] + (max-min)/n;
56       return bins;
57     }
58     Float_t dp   = n / TMath::Log10(max / min);
59     Float_t pmin = TMath::Log10(min);
60     for (Int_t i = 1; i < n+1; i++) {
61       Float_t p = pmin + i / dp;
62       bins[i]   = TMath::Power(10, p);
63     }
64     return bins;
65   }
66   //__________________________________________________________________
67   DrawHits(const char* pdgName="pi+",
68            Int_t m=1000, Double_t emin=1, Double_t emax=1000, 
69            Int_t n=900, Double_t tmin=1e-2, Double_t tmax=1e3) 
70     : fPdg(0), 
71       fRho(2.33),
72       fBetaGammaMip(3.4601)
73   { 
74     AddLoad(kKinematics);
75     AddLoad(kHits);
76     TDatabasePDG* pdgDB = TDatabasePDG::Instance();
77     fPdg                = pdgDB->GetParticle(pdgName);
78     if (!fPdg) Warning("DrawHits", "Particle %s not found", pdgName);
79
80     TArrayF tkine(MakeLogScale(n, tmin, tmax));
81     TArrayF eloss(MakeLogScale(m, emin, emax));
82     TString name("elossVsPMQ");
83     TString title(Form("#Delta E/#Delta x / q^{2} vs. p/m, %s", 
84                        (pdgName ? pdgName : "")));
85     fElossVsPMQ = new TH2D(name.Data(), title.Data(), 
86                            tkine.fN-1, tkine.fArray, 
87                            eloss.fN-1, eloss.fArray);
88     fElossVsPMQ->SetXTitle("p/(mq^{2})=#beta#gamma/q^{2}");
89     fElossVsPMQ->SetYTitle("#Delta E/#Delta x / q^{2} [MeV/cm]");
90     fEloss = new TH1D("eloss", "#Delta E/#Delta x / q^{2}", 
91                       eloss.fN-1, eloss.fArray);
92     fEloss->SetFillColor(2);
93     fEloss->SetFillStyle(3001);
94     fEloss->SetXTitle("#Delta E/#Delta x / q^{2} [MeV/cm]");
95   }
96   //__________________________________________________________________
97   Bool_t ProcessHit(AliFMDHit* hit, TParticle* p) 
98   {
99     if (!hit) {
100       std::cout << "No hit" << std::endl;
101       return kFALSE;
102     }
103
104     if (!p) {
105       std::cout << "No track" << std::endl;
106       return kFALSE;
107     }
108     if (!p->IsPrimary()) return kTRUE;
109     if (hit->IsStop()) return kTRUE;
110
111     Float_t x = hit->P();
112     Float_t y = hit->Edep()/hit->Length();
113     Float_t q = hit->Q() / 3.;
114     Float_t m = hit->M();
115     if (m == 0 || q == 0) return kTRUE;
116
117     x /= hit->M();
118     y /= q * q;
119     fElossVsPMQ->Fill(x, y);
120     fEloss->Fill(y);
121
122     return kTRUE;
123   }
124   //__________________________________________________________________
125   Bool_t Finish()
126   {
127     gStyle->SetPalette(1);
128     // gStyle->SetOptTitle(0);
129     gStyle->SetTitleBorderSize(1);
130     gStyle->SetTitleFillColor(0);
131     gStyle->SetCanvasColor(0);
132     gStyle->SetCanvasColor(0);
133     gStyle->SetCanvasBorderSize(0);
134     gStyle->SetPadColor(0);
135     gStyle->SetPadBorderSize(0);
136     TCanvas* c = new TCanvas("elossVsP", "Energy loss versus momentum");
137     c->SetLogy();
138     c->SetLogx();
139
140     TString title(fElossVsPMQ->GetTitle());
141     title.Append(Form(", %d events", fEventCount));
142     fElossVsPMQ->SetTitle(title.Data());
143     fElossVsPMQ->SetStats(kFALSE);
144     fElossVsPMQ->Draw("AXIS");
145     fElossVsPMQ->Draw("ACOLZ same");
146     TGraph*  mate    = FromGFMATE();
147     TGraph*  bb      = FromRPPFull();
148     TGraph*  nodelta = FromRPPNoDelta();
149     TGraph*  norad   = FromRPPNoRad();
150     TGraph*  mean    = FromRPPMean();
151     // fElossVsPMQ->Draw("ACOLZ same");
152     TLegend* l       = new TLegend(.5, .6, .89, .89);
153     // l->AddEntry(fElossVsPMQ, fElossVsPMQ->GetTitle(), "pf");
154     l->AddEntry(mate,        mate->GetTitle(),    "lf");
155     l->AddEntry(bb,          bb->GetTitle(),      "l");
156     l->AddEntry(nodelta,     nodelta->GetTitle(), "l");
157     l->AddEntry(norad,       norad->GetTitle(),   "l");
158     l->AddEntry(mean,        mean->GetTitle(),    "l");
159     l->Draw("same");
160     Double_t min = fElossVsPMQ->GetYaxis()->GetFirst();
161     TArrow* a = new TArrow(fBetaGammaMip, min, fBetaGammaMip, 35*min, 03, "<|");
162     a->SetAngle(30);
163     a->Draw("same");
164     TLatex* t = new TLatex(fBetaGammaMip, 40*min, "Minimum Ionising");
165     t->SetTextSize(0.04);
166     t->SetTextAlign(21);
167     t->Draw("same");
168     c->Modified();
169     c->Update();
170     c->cd();
171
172     c = new TCanvas("eloss", "Energy loss per unit material");
173     // c->SetLogx();
174     c->SetLogy();
175     fEloss->GetXaxis()->SetRangeUser(1, 10);
176     fEloss->Fit("landau", "", "", 1, 10);
177     TF1* land = fEloss->GetFunction("landau");
178     land->SetLineWidth(2);
179     Double_t max = fEloss->GetMaximum();
180     TGraph* resp  = GetResp();
181     Double_t* x   = resp->GetX();
182     Double_t* y   = resp->GetY();
183     TGraph*   g   = new TGraph(resp->GetN());
184     TGraph*   co  = GetCorr();
185     std::cout << "Correction factor: " << co->Eval(fBetaGammaMip) << std::endl;
186     Double_t  xs  = fRho; // * 1.19; // / 
187     for (Int_t i = 0; i < g->GetN(); i++) 
188       g->SetPoint(i, x[i] * xs, y[i] * max);
189     g->Draw("C same");
190
191     l = new TLegend(.6, .6, .89, .89);
192     l->AddEntry(fEloss, fEloss->GetTitle(), "lf");
193     l->AddEntry(land,   "Landau fit", "l");
194     l->AddEntry(resp,   "f(#Delta_{p}/x) [RPP fig 27.8]");
195     l->Draw("same");
196
197     return kTRUE;
198   }
199
200   /** Scale a graph by density (multiply) and mass (divide). 
201       @param graph Graph to scale 
202       @param density If @c true, scale by the Si density
203       ($\rho=2.33$/cm^3$).  The y axis is assumed to have units of
204       $MeVg^{-1}cm^2$. 
205       @param mass Mass to scale with. The x axis is assumed to be the
206       kinetic energy of the particles in units of $GeV$.  */
207   void ScaleGraph(TGraph* graph, bool density=true, double mass=1) 
208   {
209       Double_t*      x   = graph->GetX();
210       Double_t*      y   = graph->GetY();
211       const Double_t rho = (density ? fRho : 1);
212       for (Int_t i = 0; i < graph->GetN(); i++) 
213         graph->SetPoint(i, x[i] / mass, y[i] * rho); 
214   }    
215   /** Draw pure Bethe-Bloc from Review of Particle Physics, fig. 27.1 
216       @return TGraph object */ 
217   TGraph* FromRPPFull() 
218   {
219     static TGraph* graph = 0;
220     if (!graph) { 
221       graph = new TGraph(20);
222       graph->GetHistogram()->SetXTitle("#beta#gamma");
223       graph->GetHistogram()->SetYTitle("#Delta E/#Delta x [MeV/cm]");
224       graph->SetFillColor(0);
225       graph->SetLineColor(2);
226       graph->SetLineStyle(1);
227       graph->SetLineWidth(1);
228       graph->SetName("full_stop");
229       graph->SetTitle("Stopping (MeVcm^{2}/g) [RPP fig 27.1]");
230       graph->SetPoint(0,0.001461622,40.17542);
231       graph->SetPoint(1,0.003775053,91.28429);
232       graph->SetPoint(2,0.01178769,202.7359);
233       graph->SetPoint(3,0.01722915,212.1938);
234       graph->SetPoint(4,0.03162278,172.8318);
235       graph->SetPoint(5,0.06028646,91.28429);
236       graph->SetPoint(6,0.09506529,51.62633);
237       graph->SetPoint(7,0.433873,5.281682);
238       graph->SetPoint(8,1.255744,1.808947);
239       graph->SetPoint(9,2.393982,1.440177);
240       graph->SetPoint(10,3.499097,1.407715);
241       graph->SetPoint(11,10.92601,1.542122);
242       graph->SetPoint(12,60.28646,1.85066);
243       graph->SetPoint(13,236.3885,2.121938);
244       graph->SetPoint(14,468.0903,2.324538);
245       graph->SetPoint(15,1208.976,2.987085);
246       graph->SetPoint(16,6670.768,7.961412);
247       graph->SetPoint(17,23341.67,24.3298);
248       graph->SetPoint(18,110651.2,104.6651);
249       graph->SetPoint(19,264896.9,260.5203);
250       ScaleGraph(graph);
251     }
252     graph->Draw("C same");
253     return graph;
254   }
255
256   /** Draw pure Bethe-Bloc from Review of Particle Physics, fig. 27.1,
257       but without delta electrons  
258       @return TGraph object */ 
259   TGraph* FromRPPNoDelta() 
260   {
261     static TGraph* graph = 0;
262     if (!graph) { 
263       graph = new TGraph(20);
264       graph->SetName("stop_nodelta");
265       graph->SetTitle("Stopping w/o #delta's [RPP fig 27.1]");
266       graph->GetHistogram()->SetYTitle("(MeVcm^{2}/g)");
267       graph->GetHistogram()->SetXTitle("#beta#gamma");
268       graph->SetFillColor(0);
269       graph->SetLineColor(3);
270       graph->SetLineStyle(1);
271       graph->SetLineWidth(1);
272       graph->SetPoint(0,0.001461622,40.17542);
273       graph->SetPoint(1,0.003775053,91.28429);
274       graph->SetPoint(2,0.01178769,202.7359);
275       graph->SetPoint(3,0.01722915,212.1938);
276       graph->SetPoint(4,0.03162278,172.8318);
277       graph->SetPoint(5,0.06028646,91.28429);
278       graph->SetPoint(6,0.09506529,51.62633);
279       graph->SetPoint(7,0.433873,5.281682);
280       graph->SetPoint(8,1.255744,1.808947);
281       graph->SetPoint(9,2.304822,1.473387);
282       graph->SetPoint(10,3.921088,1.473387);
283       graph->SetPoint(11,8.064796,1.614064);
284       graph->SetPoint(12,26.15667,1.936996);
285       graph->SetPoint(13,264.8969,2.489084);
286       graph->SetPoint(14,544.8334,2.665278);
287       graph->SetPoint(15,1163.949,2.853945);
288       graph->SetPoint(16,5312.204,3.19853);
289       graph->SetPoint(17,15374.93,3.424944);
290       graph->SetPoint(18,49865.73,3.667384);
291       graph->SetPoint(19,634158.5,4.110185);
292       ScaleGraph(graph);
293     }
294     graph->Draw("C same");
295     return graph;
296   }
297
298   /** Draw pure Bethe-Bloc from Review of Particle Physics, fig. 27.1,
299       but without delta electrons  
300       @return TGraph object */ 
301   TGraph* FromRPPNoRad() 
302   {
303     static TGraph* graph = 0;
304     if (!graph) { 
305       graph = new TGraph(18);
306       graph->SetName("norad_stop");
307       graph->SetTitle("Stopping w/o radiative loss [RPP fig. 27.1]");
308       graph->GetHistogram()->SetYTitle("(MeVcm^{2}/g)");
309       graph->GetHistogram()->SetXTitle("#beta#gamma");
310       graph->SetFillColor(0);
311       graph->SetLineStyle(1);
312       graph->SetLineColor(4);
313       graph->SetLineWidth(1);
314       graph->SetPoint(0,0.001,24.3298);
315       graph->SetPoint(1,0.003117649,74.35105);
316       graph->SetPoint(2,0.008675042,172.8318);
317       graph->SetPoint(3,0.01782497,212.1938);
318       graph->SetPoint(4,0.02704573,189.3336);
319       graph->SetPoint(5,0.07481082,70.29816);
320       graph->SetPoint(6,0.3300035,8.524974);
321       graph->SetPoint(7,0.819559,2.489084);
322       graph->SetPoint(8,1.447084,1.651284);
323       graph->SetPoint(9,2.555097,1.440177);
324       graph->SetPoint(10,4.026598,1.407715);
325       graph->SetPoint(11,32.38084,1.728318);
326       graph->SetPoint(12,97.19733,1.893336);
327       graph->SetPoint(13,1732.539,2.170869);
328       graph->SetPoint(14,11098.58,2.324538);
329       graph->SetPoint(15,32075.46,2.378141);
330       graph->SetPoint(16,221655.8,2.546482);
331       graph->SetPoint(17,593830.6,2.605203);
332       ScaleGraph(graph);
333     }
334     graph->Draw("C same");
335     return graph;
336   }
337
338   /** Draw pure Bethe-Bloc from Review of Particle Physics, fig. 27.6 
339       @return TGraph object */ 
340   TGraph* FromRPPMean() 
341   {
342     static TGraph* graph = 0;
343     if (!graph) { 
344       graph = new TGraph(12);
345       graph->SetName("mean_eloss");
346       graph->SetTitle("Mean #Delta E/#Delta x - electronic only  [RPP fig. 27.6]");
347       graph->GetHistogram()->SetYTitle("(MeVcm^{2}/g)");
348       graph->GetHistogram()->SetXTitle("#mu E_{kin} (GeV)");
349       graph->SetFillColor(1);
350       graph->SetLineStyle(1);
351       graph->SetLineWidth(1);
352       graph->SetLineColor(1);
353       graph->SetMarkerStyle(21);
354       graph->SetMarkerSize(0.6);
355       graph->SetPoint(0,0.1,1.346561);
356       graph->SetPoint(1,0.1435819,1.230159);
357       graph->SetPoint(2,0.2061576,1.156085);
358       graph->SetPoint(3,0.3698076,1.124339);
359       graph->SetPoint(4,0.4620113,1.124339);
360       graph->SetPoint(5,0.8521452,1.145503);
361       graph->SetPoint(6,1.909707,1.177249);
362       graph->SetPoint(7,4.048096,1.198413);
363       graph->SetPoint(8,12.66832,1.219577);
364       graph->SetPoint(9,48.17031,1.230159);
365       graph->SetPoint(10,285.8863,1.230159);
366       graph->SetPoint(11,894.6674,1.230159);
367       const Double_t m   = 0.10566; // Muon 
368       ScaleGraph(graph, true, m);
369     }
370     graph->Draw("C same");
371     return graph;
372   }
373
374   /** Draw energy loss as obtained from GEANT 3.21 GFMATE. 
375       @return TGraph object */
376   TGraph* FromGFMATE() 
377   {
378     static TGraphErrors* gre = 0;
379     if (!gre) {
380       gre = new TGraphErrors(91);
381       gre->SetName("ELOSS");
382       gre->SetTitle("Energy loss 300#mu Si [GFMATE]");
383       gre->GetHistogram()->SetXTitle("#beta#gamma");
384       gre->GetHistogram()->SetYTitle("#Delta E/#Delta x [MeV/cm]");
385       gre->SetFillColor(6);
386       gre->SetFillStyle(3002);
387       gre->SetLineColor(1);
388       gre->SetLineStyle(1);
389       gre->SetLineWidth(1);
390       gre->SetPoint(0,7.16486e-05,1218.84);
391       gre->SetPoint(1,9.25378e-05,1221.38);
392       gre->SetPoint(2,0.000119517,1180.12);
393       gre->SetPoint(3,0.000154362,1100.31);
394       gre->SetPoint(4,0.000199367,996.621);
395       gre->SetPoint(5,0.000257492,886.005);
396       gre->SetPoint(6,0.000332563,780.483);
397       gre->SetPoint(7,0.000429522,684.927);
398       gre->SetPoint(8,0.000554749,599.407);
399       gre->SetPoint(9,0.000716486,522.375);
400       gre->SetPoint(10,0.000925378,452.497);
401       gre->SetPoint(11,0.00119517,389.101);
402       gre->SetPoint(12,0.00154362,331.974);
403       gre->SetPoint(13,0.00199367,280.969);
404       gre->SetPoint(14,0.00257492,235.689);
405       gre->SetPoint(15,0.00332564,196.156);
406       gre->SetPoint(16,0.00429522,162.402);
407       gre->SetPoint(17,0.00554749,133.87);
408       gre->SetPoint(18,0.00716486,109.959);
409       gre->SetPoint(19,0.00925378,90.2035);
410       gre->SetPoint(20,0.0119517,74.1317);
411       gre->SetPoint(21,0.0154362,60.8988);
412       gre->SetPoint(22,0.0199367,49.9915);
413       gre->SetPoint(23,0.0257492,40.9812);
414       gre->SetPoint(24,0.0332564,33.5739);
415       gre->SetPoint(25,0.0429522,27.5127);
416       gre->SetPoint(26,0.0554749,22.5744);
417       gre->SetPoint(27,0.0716486,18.5674);
418       gre->SetPoint(28,0.0925378,15.3292);
419       gre->SetPoint(29,0.119517,12.7231);
420       gre->SetPoint(30,0.154362,10.6352);
421       gre->SetPoint(31,0.199367,8.97115);
422       gre->SetPoint(32,0.257492,7.65358);
423       gre->SetPoint(33,0.332564,6.61909);
424       gre->SetPoint(34,0.429522,5.79837);
425       gre->SetPoint(35,0.554749,5.148);
426       gre->SetPoint(36,0.716486,4.65024);
427       gre->SetPoint(37,0.925378,4.27671);
428       gre->SetPoint(38,1.19517,3.99831);
429       gre->SetPoint(39,1.54362,3.79877);
430       gre->SetPoint(40,1.99367,3.6629);
431       gre->SetPoint(41,2.57492,3.57594);
432       gre->SetPoint(42,3.32564,3.52565);
433       gre->SetPoint(43,4.29522,3.50206);
434       gre->SetPoint(44,5.54749,3.49715);
435       gre->SetPoint(45,7.16486,3.50467);
436       gre->SetPoint(46,9.25378,3.51988);
437       gre->SetPoint(47,11.9517,3.53932);
438       gre->SetPoint(48,15.4362,3.56054);
439       gre->SetPoint(49,19.9367,3.58189);
440       gre->SetPoint(50,25.7492,3.60231);
441       gre->SetPoint(51,33.2564,3.62113);
442       gre->SetPoint(52,42.9522,3.638);
443       gre->SetPoint(53,55.4749,3.65275);
444       gre->SetPoint(54,71.6486,3.66537);
445       gre->SetPoint(55,92.5378,3.67586);
446       gre->SetPoint(56,119.517,3.68433);
447       gre->SetPoint(57,154.362,3.69105);
448       gre->SetPoint(58,199.367,3.6962);
449       gre->SetPoint(59,257.492,3.69997);
450       gre->SetPoint(60,332.564,3.70257);
451       gre->SetPoint(61,429.522,3.70421);
452       gre->SetPoint(62,554.749,3.70511);
453       gre->SetPoint(63,716.486,3.7055);
454       gre->SetPoint(64,925.378,3.70559);
455       gre->SetPoint(65,1195.17,3.70558);
456       gre->SetPoint(66,1543.62,3.70557);
457       gre->SetPoint(67,1993.67,3.70555);
458       gre->SetPoint(68,2574.92,3.70553);
459       gre->SetPoint(69,3325.64,3.70552);
460       gre->SetPoint(70,4295.22,3.7055);
461       gre->SetPoint(71,5547.49,3.70548);
462       gre->SetPoint(72,7164.86,3.70547);
463       gre->SetPoint(73,9253.78,3.70545);
464       gre->SetPoint(74,11951.7,3.70544);
465       gre->SetPoint(75,15436.2,3.70544);
466       gre->SetPoint(76,19936.7,3.70544);
467       gre->SetPoint(77,25749.2,3.70544);
468       gre->SetPoint(78,33256.4,3.70544);
469       gre->SetPoint(79,42952.2,3.70544);
470       gre->SetPoint(80,55474.9,3.70544);
471       gre->SetPoint(81,71648.6,3.70544);
472       gre->SetPoint(82,92537.8,3.70544);
473       gre->SetPoint(83,119517,3.70544);
474       gre->SetPoint(84,154362,3.70544);
475       gre->SetPoint(85,199367,3.70544);
476       gre->SetPoint(86,257492,3.70544);
477       gre->SetPoint(87,332563,3.70544);
478       gre->SetPoint(88,429522,3.70544);
479       gre->SetPoint(89,554749,3.70544);
480       gre->SetPoint(90,716486,3.70544);
481       // Double_t* x = gre->GetX();
482       Double_t* y = gre->GetY();
483       for (Int_t i = 0; i < gre->GetN(); i++) 
484         gre->SetPointError(i, 0, 2 * 0.1 * y[i]); // ! 1 sigma
485     }
486     gre->DrawClone("c3 same");
487     return gre;
488   }
489
490   /** Get the response functin @f$ f(\Delta_p/x)@f$ from Review of
491       Particle Physics (fig. 27.2).  It is scaled to the value at
492       MPV. */ 
493   TGraph* GetResp()
494   {
495     static TGraph*  graph = 0;
496     if (!graph) {
497       graph = new TGraph(16);
498       graph->SetName("si_resp");
499       graph->SetTitle("f(#Delta_{p}/x) scaled to the MPV value ");
500       graph->GetHistogram()->SetXTitle("#Delta_{p}/x (MeVcm^{2}/g)");
501       graph->GetHistogram()->SetYTitle("f(#Delta_{p}/x)");
502       graph->SetFillColor(1);
503       graph->SetMarkerStyle(21);
504       graph->SetMarkerSize(0.6);
505       graph->SetPoint(0,0.8115124,0.009771987);
506       graph->SetPoint(1,0.9198646,0.228013);
507       graph->SetPoint(2,0.996614,0.5895765);
508       graph->SetPoint(3,1.041761,0.8241042);
509       graph->SetPoint(4,1.059819,0.8794788);
510       graph->SetPoint(5,1.077878,0.9348534);
511       graph->SetPoint(6,1.100451,0.980456);
512       graph->SetPoint(7,1.141084,0.9967427);
513       graph->SetPoint(8,1.204289,0.9153094);
514       graph->SetPoint(9,1.276524,0.742671);
515       graph->SetPoint(10,1.402935,0.465798);
516       graph->SetPoint(11,1.515801,0.3029316);
517       graph->SetPoint(12,1.73702,0.1465798);
518       graph->SetPoint(13,1.985327,0.08143322);
519       graph->SetPoint(14,2.301354,0.04234528);
520       graph->SetPoint(15,2.56772,0.02931596);
521     }
522     return graph;
523   }
524
525   /** Get the correction to Bethe-Bloc from Review of Particle Physics
526       (fig 27.8). 
527   */
528   TGraph* GetCorr() 
529   {
530     static TGraph* graph = 0;
531     if (!graph) {
532       graph = new TGraph(14);
533       graph->SetName("graph");
534       graph->SetTitle("(#Delta_{p}/x)/(dE/dx)|_{mip} for 320#mu Si");
535       graph->GetHistogram()->SetXTitle("#beta#gamma = p/m");
536       graph->SetFillColor(1);
537       graph->SetLineColor(7);
538       graph->SetMarkerStyle(21);
539       graph->SetMarkerSize(0.6);
540       graph->SetPoint(0,1.196058,0.9944915);
541       graph->SetPoint(1,1.28502,0.9411017);
542       graph->SetPoint(2,1.484334,0.8559322);
543       graph->SetPoint(3,1.984617,0.7491525);
544       graph->SetPoint(4,2.658367,0.6983051);
545       graph->SetPoint(5,3.780227,0.6779661);
546       graph->SetPoint(6,4.997358,0.6741525);
547       graph->SetPoint(7,8.611026,0.684322);
548       graph->SetPoint(8,15.28296,0.6995763);
549       graph->SetPoint(9,41.54516,0.7186441);
550       graph->SetPoint(10,98.91461,0.7288136);
551       graph->SetPoint(11,203.2734,0.7326271);
552       graph->SetPoint(12,505.6421,0.7338983);
553       graph->SetPoint(13,896.973,0.7338983);
554     }
555     return graph;
556   }
557   
558   ClassDef(DrawHits,0);
559 };
560
561 //____________________________________________________________________
562 //
563 // EOF
564 //