Script to get # of dead channels from OCDB
[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 <AliStack.h>
18 #include <iostream>
19 #include <TStyle.h>
20 #include <TArrayF.h>
21 #include <TParticle.h>
22 #include <TCanvas.h>
23 #include <TGraphErrors.h>
24 #include <TDatabasePDG.h>
25 #include <TParticlePDG.h>
26 #include <TLegend.h>
27 #include <TArrow.h>
28 #include <TLatex.h>
29 #include <TF1.h>
30
31 /** @class DrawHits
32     @brief Draw hit energy loss
33     @code 
34     Root> .L Compile.C
35     Root> Compile("DrawHits.C")
36     Root> DrawHits c
37     Root> c.Run();
38     @endcode
39     @ingroup FMD_script
40  */
41 class DrawHits : public AliFMDInput
42 {
43 private:
44   TH2D* fElossVsPMQ; // Histogram 
45   TH1D* fEloss;
46   TH1D* fBetaGamma;
47   TParticlePDG* fPdg;
48   const Double_t fRho;
49   const Double_t fBetaGammaMip;
50 public:
51   //__________________________________________________________________
52   DrawHits(const char* pdgName="pi+",
53            Int_t m=500, Double_t emin=1, Double_t emax=1000, 
54            Int_t n=900, Double_t tmin=1e-2, Double_t tmax=1e3) 
55     : AliFMDInput("galice.root"),
56       fElossVsPMQ(0),
57       fEloss(0),
58       fBetaGamma(0),
59       fPdg(0), 
60       fRho(2.33), // 2.33),
61       fBetaGammaMip(3.4601)
62   { 
63     AddLoad(kKinematics);
64     AddLoad(kHits);
65     TDatabasePDG* pdgDB = TDatabasePDG::Instance();
66     fPdg                = pdgDB->GetParticle(pdgName);
67     if (!fPdg) Warning("DrawHits", "Particle %s not found", pdgName);
68
69     TArrayF tkine(MakeLogScale(n,    tmin, tmax));
70     TArrayF betag(MakeLogScale(n/10, tmin, tmax));
71     TArrayF eloss(MakeLogScale(m,    emin, emax));
72     TString name("elossVsPMQ");
73     TString title(Form("#Delta E/#Delta x / q^{2} vs. p/m, %s", 
74                        (pdgName ? pdgName : "")));
75     fElossVsPMQ = new TH2D(name.Data(), title.Data(), 
76                            tkine.fN-1, tkine.fArray, 
77                            eloss.fN-1, eloss.fArray);
78     fElossVsPMQ->SetXTitle("p/(mq^{2})=#beta#gamma/q^{2}");
79     fElossVsPMQ->SetYTitle("#Delta E/#Delta x / q^{2} [MeV/cm]");
80     fElossVsPMQ->Sumw2();
81     fEloss = new TH1D("eloss", "#Delta E/#Delta x / q^{2}", 
82                       eloss.fN-1, eloss.fArray);
83     fEloss->SetFillColor(kRed);
84     fEloss->SetFillStyle(3001);
85     fEloss->SetXTitle("#Delta E/#Delta x / q^{2} [MeV/cm]");
86     fEloss->Sumw2();
87
88     fBetaGamma = new TH1D("betaGamma", "#beta#gamma of particles", 
89                           betag.fN-1, betag.fArray);
90     fBetaGamma->SetXTitle("#beta#gamma");
91     fBetaGamma->SetFillColor(kBlue+1);
92     fBetaGamma->SetFillStyle(3001);
93   }
94   //__________________________________________________________________
95   Bool_t ProcessHit(AliFMDHit* hit, TParticle* p) 
96   {
97     if (!hit) {
98       std::cout << "No hit" << std::endl;
99       return kFALSE;
100     }
101
102     if (!p) {
103       std::cout << "No track" << std::endl;
104       return kFALSE;
105     }
106     // if (!p->IsPrimary()) return kTRUE;
107     if (hit->IsStop()) return kTRUE;
108     if (hit->Length() == 0) { 
109       Warning("ProcessHit", "Hit in %s has 0 length", hit->GetName());
110       return kTRUE;
111     }
112     
113     Float_t q = hit->Q() / 3.;
114     Float_t m = hit->M();
115     if (m == 0 || q == 0) return kTRUE;
116
117     TLorentzVector pp;
118     p->Momentum(pp);
119     Double_t betagamma = 0;
120     Info("ProcessHit", "%s (%s) beta=%f", p->GetPDG()->GetName(), 
121          fStack->IsPhysicalPrimary(hit->Track()) ? "primary" : "secondary", 
122          pp.Beta());
123     if (pp.Beta() <= 1 && pp.Beta() >= 0) 
124       betagamma = pp.Beta() * pp.Gamma();
125     fBetaGamma->Fill(betagamma);
126 #if 0
127     if (betagamma < 10) { 
128       Info("ProcessHit", "%s (%s) beta=%f gamma=%f beta*gamma=%f", 
129            p->GetPDG()->GetName(), 
130            fStack->IsPhysicalPrimary(hit->Track()) ? "primary" : 
131            "secondary", 
132            pp.Beta(), pp.Gamma(), betagamma);
133       return kTRUE;
134     }
135 #endif
136
137     Float_t x = hit->P();
138     Float_t y = hit->Edep()/hit->Length();
139
140     x /= hit->M();
141     // y /= q * q;
142     fElossVsPMQ->Fill(x, y);
143     fEloss->Fill(y);
144     // fNHits++;
145     return kTRUE;
146   }
147   //__________________________________________________________________
148   void ShowFit(Double_t x1, Double_t y1, const char* title, 
149                TF1* f, Double_t dx=0, Double_t dy=0.05)
150   {
151     Double_t x = x1, y = y1;
152     TLatex* latex = new TLatex(x, y, title);
153     latex->SetTextFont(132);
154     latex->SetTextSize(0.8*dy);
155     latex->SetNDC();
156     latex->Draw();
157     x -= dx;
158     y -= dy;
159     const Double_t eqDx=0.1;
160     Double_t chi2 = f->GetChisquare();
161     Int_t    ndf  = f->GetNDF();
162     Double_t prob = f->GetProb();
163     latex->DrawLatex(x, y, "#chi^{2}/NDF");
164     latex->DrawLatex(x+eqDx, y, Form("= %7.4f/%3d=%5.2f (%3d%%)", 
165                                      chi2, ndf, chi2/ndf, int(100*prob)));
166     Int_t     n = f->GetNpar();
167     Double_t* p = f->GetParameters();
168     Double_t* e = f->GetParErrors();
169     for (int i = 0; i < n; i++) { 
170       x -= dx;
171       y -= dy;
172       latex->DrawLatex(x, y, f->GetParName(i));
173       latex->DrawLatex(x+eqDx, y, Form("= %7.4f", p[i]));
174       latex->DrawLatex(x+2*eqDx, y, Form("#pm %7.4f", e[i]));
175     }
176   }      
177   //__________________________________________________________________
178   Bool_t Finish()
179   {
180     gStyle->SetPalette(1);
181     // gStyle->SetOptTitle(0);
182     gStyle->SetTitleBorderSize(1);
183     gStyle->SetTitleFillColor(0);
184     gStyle->SetCanvasColor(0);
185     gStyle->SetCanvasColor(0);
186     gStyle->SetCanvasBorderSize(0);
187     gStyle->SetPadColor(0);
188     gStyle->SetPadBorderSize(0);
189     gStyle->SetOptStat(0);
190     TCanvas* c = new TCanvas("elossVsP", "Energy loss versus momentum", 
191                              1200, 800);
192     c->SetLogy();
193     c->SetLogx();
194
195     TString title(Form("%s, %d events", fElossVsPMQ->GetTitle(), fEventCount));
196     fElossVsPMQ->SetTitle(title.Data());
197     fElossVsPMQ->SetStats(kFALSE);
198     fElossVsPMQ->Draw("AXIS");
199     fElossVsPMQ->Draw("ACOLZ same");
200     TGraph*  mate    = FromGFMATE();
201     TGraph*  bb      = FromRPPFull();
202     TGraph*  nodelta = FromRPPNoDelta();
203     TGraph*  norad   = FromRPPNoRad();
204     TGraph*  mean    = FromRPPMean();
205     // fElossVsPMQ->Draw("ACOLZ same");
206     TLegend* l       = new TLegend(.5, .6, .89, .89);
207     // l->AddEntry(fElossVsPMQ, fElossVsPMQ->GetTitle(), "pf");
208     l->SetFillColor(0);
209     l->AddEntry(mate,        mate->GetTitle(),    "lf");
210     l->AddEntry(bb,          bb->GetTitle(),      "l");
211     l->AddEntry(nodelta,     nodelta->GetTitle(), "l");
212     l->AddEntry(norad,       norad->GetTitle(),   "l");
213     l->AddEntry(mean,        mean->GetTitle(),    "l");
214     l->Draw("same");
215     Double_t min = fElossVsPMQ->GetYaxis()->GetFirst();
216     TArrow* a = new TArrow(fBetaGammaMip, min, fBetaGammaMip, 35*min, 03, "<|");
217     a->SetAngle(30);
218     a->Draw("same");
219     TLatex* t = new TLatex(fBetaGammaMip, 40*min, "Minimum Ionising");
220     t->SetTextSize(0.04);
221     t->SetTextAlign(21);
222     t->Draw("same");
223     c->Modified();
224     c->Update();
225     c->cd();
226     c->SaveAs("eloss_bethe.png");
227
228     c = new TCanvas("cEloss", "Energy loss per unit material",
229                     1200, 800);
230     c->SetRightMargin(0.05);
231     c->SetTopMargin(0.05);
232     c->SetLeftMargin(0.05);
233     fEloss->SetStats(kFALSE);
234     // c->SetLogx();
235     TF1* land     = new TF1("land", "landau");
236     land->SetParNames("A", "MPV", "width");
237     land->SetLineWidth(2);
238     land->SetLineColor(kGreen+1);
239
240     TF1* landgaus = new TF1("landgaus", "gaus(0)+landau(3)");
241     landgaus->SetParNames("A", "#mu", "#sigma", "B", "MPV", "width");
242     landgaus->SetLineWidth(2);
243     landgaus->SetLineColor(kMagenta+1);
244     TGraph*  corr  = GetCorr();
245     TGraph*  resp  = GetResp();
246     if (fEloss->GetEntries() != 0) { 
247       c->SetLogy();
248       fEloss->Scale(1. / fEloss->GetEntries());
249       fEloss->GetXaxis()->SetRangeUser(1, 10);
250       fEloss->Fit(land, "+", "", 2, 10);
251       landgaus->SetParameters(land->GetParameter(0) / 100, 
252                               land->GetParameter(1), 
253                               land->GetParameter(2), 
254                               land->GetParameter(0),
255                               land->GetParameter(1), 
256                               land->GetParameter(2));
257       fEloss->Fit(landgaus, "+", "", 1, 10);
258       fEloss->DrawCopy("HIST SAME");
259     }
260     
261     // fEloss->DrawCopy("E SAME");
262     // land->Draw("same");
263     // landgaus->Draw("same");
264     Double_t max = fEloss->GetMaximum();
265     Double_t* x   = resp->GetX();
266     Double_t* y   = resp->GetY();
267     TGraph*   g   = new TGraph(resp->GetN());
268     g->SetName(Form("%sCorr", resp->GetName()));
269     g->SetTitle(resp->GetTitle());
270     g->SetLineStyle(resp->GetLineStyle());
271     g->SetLineColor(resp->GetLineColor());
272     g->SetLineWidth(resp->GetLineWidth());
273     Double_t  xs2 = corr->Eval(fBetaGammaMip);
274     Double_t xss   = 1.1;
275     Double_t  xs  = fRho * xss;
276     std::cout << "Correction factor: " << xs2 << std::endl;
277     for (Int_t i = 0; i < g->GetN(); i++) 
278       g->SetPoint(i, x[i] * xs, y[i] * max);
279     g->Draw("C same");
280     
281     l = new TLegend(.05, .6, .4, .95);
282     l->SetFillColor(0);
283     l->SetBorderSize(1);
284     l->AddEntry(fEloss, fEloss->GetTitle(), "lf");
285     if (land) 
286       l->AddEntry(land,   Form("Landau fit\t- #chi^{2}/NDF=%7.5f", 
287                                land->GetChisquare()/land->GetNDF()), "l");
288     if (landgaus) 
289       l->AddEntry(landgaus,   
290                   Form("Landau+Gauss fit\t- #chi^{2}/NDF=%7.5f", 
291                        landgaus->GetChisquare()/landgaus->GetNDF()), "l");
292     l->AddEntry(resp,   Form("f(%s#Delta/x) 320#mum Si [RPP fig 27.7]",
293                              xss != 1 ? Form("%4.2f#times", xss) : ""), 
294                 "l");
295     l->Draw("same");
296     
297     fEloss->GetYaxis()->SetRangeUser(1e-4, 100);
298     ShowFit(0.45,.9, "Landau+Gaus", landgaus, 0, 0.02);
299     ShowFit(0.70,.9, "Landau", land, 0, 0.02);
300
301     c->Modified();
302     c->Update();
303     c->cd();
304     c->SaveAs("eloss_landau.png");
305
306
307     c = new TCanvas("cBetaGamma", "beta gamma", 1200, 800);
308     c->SetLogx();
309     fBetaGamma->Draw();
310     Int_t mipbin = fBetaGamma->FindBin(fBetaGammaMip) + 1;
311     Int_t maxbin = fBetaGamma->GetNbinsX();
312     Int_t total  = fBetaGamma->Integral();
313     Int_t over   = fBetaGamma->Integral(mipbin,maxbin);
314     TH1*  res    = (TH1*)fBetaGamma->Clone("overMip");
315     res->SetFillColor(kRed+1);
316     for (int i = 0; i < mipbin; i++) 
317       res->SetBinContent(i, 0);
318     res->Draw("same");
319     std::cout << "Percentage over MIP : " << float(over) / total << std::endl;
320
321     Double_t yy = fBetaGamma->GetBinContent(mipbin) * 1.1; 
322     a = new TArrow(fBetaGammaMip, 0, fBetaGammaMip, yy, 3, "<|");
323     a->SetAngle(30);
324     a->Draw("same");
325     t = new TLatex(fBetaGammaMip, yy, "Minimum Ionising");
326     t->SetTextSize(0.04);
327     t->SetTextAlign(21);
328     t->Draw("same");
329     c->Modified();
330     c->Update();
331     c->cd();
332     c->SaveAs("eloss_betagamma.png");
333     
334     return kTRUE;
335   }
336
337   /** Scale a graph by density (multiply) and mass (divide). 
338       @param graph Graph to scale 
339       @param density If @c true, scale by the Si density
340       ($\rho=2.33$/cm^3$).  The y axis is assumed to have units of
341       $MeVg^{-1}cm^2$. 
342       @param mass Mass to scale with. The x axis is assumed to be the
343       kinetic energy of the particles in units of $GeV$.  */
344   void ScaleGraph(TGraph* graph, bool density=true, double mass=1) 
345   {
346       Double_t*      x   = graph->GetX();
347       Double_t*      y   = graph->GetY();
348       const Double_t rho = (density ? fRho : 1);
349       for (Int_t i = 0; i < graph->GetN(); i++) 
350         graph->SetPoint(i, x[i] / mass, y[i] * rho); 
351   }    
352   /** Draw pure Bethe-Bloc from Review of Particle Physics, fig. 27.1 
353       @return TGraph object */ 
354   TGraph* FromRPPFull() 
355   {
356     static TGraph* graph = 0;
357     if (!graph) { 
358       graph = new TGraph(20);
359       graph->GetHistogram()->SetXTitle("#beta#gamma");
360       graph->GetHistogram()->SetYTitle("#Delta E/#Delta x [MeV/cm]");
361       graph->SetFillColor(0);
362       graph->SetLineColor(kRed+1);
363       graph->SetLineStyle(2);
364       graph->SetLineWidth(2);
365       graph->SetName("full_stop");
366       graph->SetTitle("Stopping (MeVcm^{2}/g) [RPP fig 27.1]");
367       graph->SetPoint(0,0.001461622,40.17542);
368       graph->SetPoint(1,0.003775053,91.28429);
369       graph->SetPoint(2,0.01178769,202.7359);
370       graph->SetPoint(3,0.01722915,212.1938);
371       graph->SetPoint(4,0.03162278,172.8318);
372       graph->SetPoint(5,0.06028646,91.28429);
373       graph->SetPoint(6,0.09506529,51.62633);
374       graph->SetPoint(7,0.433873,5.281682);
375       graph->SetPoint(8,1.255744,1.808947);
376       graph->SetPoint(9,2.393982,1.440177);
377       graph->SetPoint(10,3.499097,1.407715);
378       graph->SetPoint(11,10.92601,1.542122);
379       graph->SetPoint(12,60.28646,1.85066);
380       graph->SetPoint(13,236.3885,2.121938);
381       graph->SetPoint(14,468.0903,2.324538);
382       graph->SetPoint(15,1208.976,2.987085);
383       graph->SetPoint(16,6670.768,7.961412);
384       graph->SetPoint(17,23341.67,24.3298);
385       graph->SetPoint(18,110651.2,104.6651);
386       graph->SetPoint(19,264896.9,260.5203);
387       ScaleGraph(graph);
388     }
389     graph->Draw("C same");
390     return graph;
391   }
392
393   /** Draw pure Bethe-Bloc from Review of Particle Physics, fig. 27.1,
394       but without delta electrons  
395       @return TGraph object */ 
396   TGraph* FromRPPNoDelta() 
397   {
398     static TGraph* graph = 0;
399     if (!graph) { 
400       graph = new TGraph(20);
401       graph->SetName("stop_nodelta");
402       graph->SetTitle("Stopping w/o #delta's [RPP fig 27.1]");
403       graph->GetHistogram()->SetYTitle("(MeVcm^{2}/g)");
404       graph->GetHistogram()->SetXTitle("#beta#gamma");
405       graph->SetFillColor(0);
406       graph->SetLineColor(kGreen+1);
407       graph->SetLineStyle(2);
408       graph->SetLineWidth(2);
409       graph->SetPoint(0,0.001461622,40.17542);
410       graph->SetPoint(1,0.003775053,91.28429);
411       graph->SetPoint(2,0.01178769,202.7359);
412       graph->SetPoint(3,0.01722915,212.1938);
413       graph->SetPoint(4,0.03162278,172.8318);
414       graph->SetPoint(5,0.06028646,91.28429);
415       graph->SetPoint(6,0.09506529,51.62633);
416       graph->SetPoint(7,0.433873,5.281682);
417       graph->SetPoint(8,1.255744,1.808947);
418       graph->SetPoint(9,2.304822,1.473387);
419       graph->SetPoint(10,3.921088,1.473387);
420       graph->SetPoint(11,8.064796,1.614064);
421       graph->SetPoint(12,26.15667,1.936996);
422       graph->SetPoint(13,264.8969,2.489084);
423       graph->SetPoint(14,544.8334,2.665278);
424       graph->SetPoint(15,1163.949,2.853945);
425       graph->SetPoint(16,5312.204,3.19853);
426       graph->SetPoint(17,15374.93,3.424944);
427       graph->SetPoint(18,49865.73,3.667384);
428       graph->SetPoint(19,634158.5,4.110185);
429       ScaleGraph(graph);
430     }
431     graph->Draw("C same");
432     return graph;
433   }
434
435   /** Draw pure Bethe-Bloc from Review of Particle Physics, fig. 27.1,
436       but without delta electrons  
437       @return TGraph object */ 
438   TGraph* FromRPPNoRad() 
439   {
440     static TGraph* graph = 0;
441     if (!graph) { 
442       graph = new TGraph(18);
443       graph->SetName("norad_stop");
444       graph->SetTitle("Stopping w/o radiative loss [RPP fig. 27.1]");
445       graph->GetHistogram()->SetYTitle("(MeVcm^{2}/g)");
446       graph->GetHistogram()->SetXTitle("#beta#gamma");
447       graph->SetFillColor(0);
448       graph->SetLineColor(kBlue+1);
449       graph->SetLineWidth(2);
450       graph->SetLineStyle(2);
451       graph->SetPoint(0,0.001,24.3298);
452       graph->SetPoint(1,0.003117649,74.35105);
453       graph->SetPoint(2,0.008675042,172.8318);
454       graph->SetPoint(3,0.01782497,212.1938);
455       graph->SetPoint(4,0.02704573,189.3336);
456       graph->SetPoint(5,0.07481082,70.29816);
457       graph->SetPoint(6,0.3300035,8.524974);
458       graph->SetPoint(7,0.819559,2.489084);
459       graph->SetPoint(8,1.447084,1.651284);
460       graph->SetPoint(9,2.555097,1.440177);
461       graph->SetPoint(10,4.026598,1.407715);
462       graph->SetPoint(11,32.38084,1.728318);
463       graph->SetPoint(12,97.19733,1.893336);
464       graph->SetPoint(13,1732.539,2.170869);
465       graph->SetPoint(14,11098.58,2.324538);
466       graph->SetPoint(15,32075.46,2.378141);
467       graph->SetPoint(16,221655.8,2.546482);
468       graph->SetPoint(17,593830.6,2.605203);
469       ScaleGraph(graph);
470     }
471     graph->Draw("C same");
472     return graph;
473   }
474
475   /** Draw pure Bethe-Bloc from Review of Particle Physics, fig. 27.6 
476       @return TGraph object */ 
477   TGraph* FromRPPMean() 
478   {
479     static TGraph* graph = 0;
480     if (!graph) { 
481       graph = new TGraph(12);
482       graph->SetName("mean_eloss");
483       graph->SetTitle("Mean #Delta E/#Delta x - "
484                       "electronic only  [RPP fig. 27.6]");
485       graph->GetHistogram()->SetYTitle("(MeVcm^{2}/g)");
486       graph->GetHistogram()->SetXTitle("#mu E_{kin} (GeV)");
487       graph->SetFillColor(0);
488       graph->SetLineStyle(2);
489       graph->SetLineWidth(2);
490       graph->SetLineColor(kMagenta+1);
491       graph->SetMarkerStyle(21);
492       graph->SetMarkerSize(0.6);
493       graph->SetPoint(0,0.1,1.346561);
494       graph->SetPoint(1,0.1435819,1.230159);
495       graph->SetPoint(2,0.2061576,1.156085);
496       graph->SetPoint(3,0.3698076,1.124339);
497       graph->SetPoint(4,0.4620113,1.124339);
498       graph->SetPoint(5,0.8521452,1.145503);
499       graph->SetPoint(6,1.909707,1.177249);
500       graph->SetPoint(7,4.048096,1.198413);
501       graph->SetPoint(8,12.66832,1.219577);
502       graph->SetPoint(9,48.17031,1.230159);
503       graph->SetPoint(10,285.8863,1.230159);
504       graph->SetPoint(11,894.6674,1.230159);
505       const Double_t m   = 0.10566; // Muon 
506       ScaleGraph(graph, true, m);
507     }
508     graph->Draw("C same");
509     return graph;
510   }
511
512   /** Draw energy loss as obtained from GEANT 3.21 GFMATE. 
513       @return TGraph object */
514   TGraph* FromGFMATE() 
515   {
516     static TGraphErrors* gre = 0;
517     if (!gre) {
518       gre = new TGraphErrors(91);
519       gre->SetName("ELOSS");
520       gre->SetTitle("Energy loss 300#mu Si [GFMATE]");
521       gre->GetHistogram()->SetXTitle("#beta#gamma");
522       gre->GetHistogram()->SetYTitle("#Delta E/#Delta x [MeV/cm]");
523       gre->SetFillColor(kGray);
524       gre->SetFillStyle(3001); // 0 /* 3002 */);
525       gre->SetLineColor(kGray+1);
526       gre->SetLineStyle(1);
527       gre->SetLineWidth(2);
528       gre->SetPoint(0,7.16486e-05,1218.84);
529       gre->SetPoint(1,9.25378e-05,1221.38);
530       gre->SetPoint(2,0.000119517,1180.12);
531       gre->SetPoint(3,0.000154362,1100.31);
532       gre->SetPoint(4,0.000199367,996.621);
533       gre->SetPoint(5,0.000257492,886.005);
534       gre->SetPoint(6,0.000332563,780.483);
535       gre->SetPoint(7,0.000429522,684.927);
536       gre->SetPoint(8,0.000554749,599.407);
537       gre->SetPoint(9,0.000716486,522.375);
538       gre->SetPoint(10,0.000925378,452.497);
539       gre->SetPoint(11,0.00119517,389.101);
540       gre->SetPoint(12,0.00154362,331.974);
541       gre->SetPoint(13,0.00199367,280.969);
542       gre->SetPoint(14,0.00257492,235.689);
543       gre->SetPoint(15,0.00332564,196.156);
544       gre->SetPoint(16,0.00429522,162.402);
545       gre->SetPoint(17,0.00554749,133.87);
546       gre->SetPoint(18,0.00716486,109.959);
547       gre->SetPoint(19,0.00925378,90.2035);
548       gre->SetPoint(20,0.0119517,74.1317);
549       gre->SetPoint(21,0.0154362,60.8988);
550       gre->SetPoint(22,0.0199367,49.9915);
551       gre->SetPoint(23,0.0257492,40.9812);
552       gre->SetPoint(24,0.0332564,33.5739);
553       gre->SetPoint(25,0.0429522,27.5127);
554       gre->SetPoint(26,0.0554749,22.5744);
555       gre->SetPoint(27,0.0716486,18.5674);
556       gre->SetPoint(28,0.0925378,15.3292);
557       gre->SetPoint(29,0.119517,12.7231);
558       gre->SetPoint(30,0.154362,10.6352);
559       gre->SetPoint(31,0.199367,8.97115);
560       gre->SetPoint(32,0.257492,7.65358);
561       gre->SetPoint(33,0.332564,6.61909);
562       gre->SetPoint(34,0.429522,5.79837);
563       gre->SetPoint(35,0.554749,5.148);
564       gre->SetPoint(36,0.716486,4.65024);
565       gre->SetPoint(37,0.925378,4.27671);
566       gre->SetPoint(38,1.19517,3.99831);
567       gre->SetPoint(39,1.54362,3.79877);
568       gre->SetPoint(40,1.99367,3.6629);
569       gre->SetPoint(41,2.57492,3.57594);
570       gre->SetPoint(42,3.32564,3.52565);
571       gre->SetPoint(43,4.29522,3.50206);
572       gre->SetPoint(44,5.54749,3.49715);
573       gre->SetPoint(45,7.16486,3.50467);
574       gre->SetPoint(46,9.25378,3.51988);
575       gre->SetPoint(47,11.9517,3.53932);
576       gre->SetPoint(48,15.4362,3.56054);
577       gre->SetPoint(49,19.9367,3.58189);
578       gre->SetPoint(50,25.7492,3.60231);
579       gre->SetPoint(51,33.2564,3.62113);
580       gre->SetPoint(52,42.9522,3.638);
581       gre->SetPoint(53,55.4749,3.65275);
582       gre->SetPoint(54,71.6486,3.66537);
583       gre->SetPoint(55,92.5378,3.67586);
584       gre->SetPoint(56,119.517,3.68433);
585       gre->SetPoint(57,154.362,3.69105);
586       gre->SetPoint(58,199.367,3.6962);
587       gre->SetPoint(59,257.492,3.69997);
588       gre->SetPoint(60,332.564,3.70257);
589       gre->SetPoint(61,429.522,3.70421);
590       gre->SetPoint(62,554.749,3.70511);
591       gre->SetPoint(63,716.486,3.7055);
592       gre->SetPoint(64,925.378,3.70559);
593       gre->SetPoint(65,1195.17,3.70558);
594       gre->SetPoint(66,1543.62,3.70557);
595       gre->SetPoint(67,1993.67,3.70555);
596       gre->SetPoint(68,2574.92,3.70553);
597       gre->SetPoint(69,3325.64,3.70552);
598       gre->SetPoint(70,4295.22,3.7055);
599       gre->SetPoint(71,5547.49,3.70548);
600       gre->SetPoint(72,7164.86,3.70547);
601       gre->SetPoint(73,9253.78,3.70545);
602       gre->SetPoint(74,11951.7,3.70544);
603       gre->SetPoint(75,15436.2,3.70544);
604       gre->SetPoint(76,19936.7,3.70544);
605       gre->SetPoint(77,25749.2,3.70544);
606       gre->SetPoint(78,33256.4,3.70544);
607       gre->SetPoint(79,42952.2,3.70544);
608       gre->SetPoint(80,55474.9,3.70544);
609       gre->SetPoint(81,71648.6,3.70544);
610       gre->SetPoint(82,92537.8,3.70544);
611       gre->SetPoint(83,119517,3.70544);
612       gre->SetPoint(84,154362,3.70544);
613       gre->SetPoint(85,199367,3.70544);
614       gre->SetPoint(86,257492,3.70544);
615       gre->SetPoint(87,332563,3.70544);
616       gre->SetPoint(88,429522,3.70544);
617       gre->SetPoint(89,554749,3.70544);
618       gre->SetPoint(90,716486,3.70544);
619       // Double_t* x = gre->GetX();
620       Double_t* y = gre->GetY();
621       for (Int_t i = 0; i < gre->GetN(); i++) 
622         gre->SetPointError(i, 0, 2 * 0.1 * y[i]); // ! 1 sigma
623     }
624     gre->Draw("c3 same");
625     return gre;
626   }
627
628   /** Get the response functin @f$ f(\Delta_p/x)@f$ from Review of
629       Particle Physics (fig. 27.7).  It is scaled to the value at
630       MPV. */ 
631   TGraph* GetResp()
632   {
633     static TGraph*  graph = 0;
634     if (!graph) {
635       graph = new TGraph;
636       graph->SetName("si_resp");
637       graph->SetTitle("f(#Delta/x) scaled to the MPV value ");
638       graph->GetHistogram()->SetXTitle("#Delta/x (MeVcm^{2}/g)");
639       graph->GetHistogram()->SetYTitle("f(#Delta/x)");
640       graph->SetLineColor(kBlue+1);
641       graph->SetLineWidth(2);
642       graph->SetFillColor(kBlue+1);
643       graph->SetMarkerStyle(21);
644       graph->SetMarkerSize(0.6);
645 #if 0
646       // Figure 27.7 or Review of Particle physics - Straggeling function in 
647       // silicon of 500 MeV Pions, normalised to unity at the most probable 
648       // value.   
649       graph->SetPoint(0,0.808094,0.00377358);
650       graph->SetPoint(1,0.860313,0.0566038);
651       graph->SetPoint(2,0.891645,0.116981);
652       graph->SetPoint(3,0.912533,0.181132);
653       graph->SetPoint(4,0.928198,0.260377);
654       graph->SetPoint(5,0.938642,0.320755);
655       graph->SetPoint(6,0.954308,0.377358);
656       graph->SetPoint(7,0.964752,0.433962);
657       graph->SetPoint(8,0.975196,0.490566);
658       graph->SetPoint(9,0.98564,0.550943);
659       graph->SetPoint(10,0.996084,0.611321);
660       graph->SetPoint(11,1.00653,0.667925);
661       graph->SetPoint(12,1.02219,0.732075);
662       graph->SetPoint(13,1.03264,0.784906);
663       graph->SetPoint(14,1.0483,0.845283);
664       graph->SetPoint(15,1.06397,0.901887);
665       graph->SetPoint(16,1.09008,0.958491);
666       graph->SetPoint(17,1.10574,0.984906);
667       graph->SetPoint(18,1.13708,1);
668       graph->SetPoint(19,1.13708,1);
669       graph->SetPoint(20,1.15796,0.988679);
670       graph->SetPoint(21,1.17363,0.966038);
671       graph->SetPoint(22,1.19974,0.916981);
672       graph->SetPoint(23,1.2154,0.89434);
673       graph->SetPoint(24,1.23629,0.837736);
674       graph->SetPoint(25,1.2624,0.784906);
675       graph->SetPoint(26,1.28329,0.724528);
676       graph->SetPoint(27,1.3094,0.664151);
677       graph->SetPoint(28,1.32507,0.611321);
678       graph->SetPoint(29,1.3564,0.550943);
679       graph->SetPoint(30,1.41384,0.445283);
680       graph->SetPoint(31,1.44517,0.392453);
681       graph->SetPoint(32,1.48695,0.335849);
682       graph->SetPoint(33,1.52872,0.286792);
683       graph->SetPoint(34,1.58094,0.237736);
684       graph->SetPoint(35,1.63838,0.196226);
685       graph->SetPoint(36,1.68016,0.169811);
686       graph->SetPoint(37,1.75326,0.135849);
687       graph->SetPoint(38,1.81593,0.113208);
688       graph->SetPoint(39,1.89426,0.0981132);
689       graph->SetPoint(40,1.96214,0.0830189);
690       graph->SetPoint(41,2.0718,0.0641509);
691       graph->SetPoint(42,2.19191,0.0490566);
692       graph->SetPoint(43,2.31723,0.0415094);
693       graph->SetPoint(44,2.453,0.0301887);
694       graph->SetPoint(45,2.53133,0.0264151);
695       graph->SetPoint(46,2.57833,0.0264151);
696 #else
697       graph->SetPoint(0,0.8115124,0.009771987);
698       graph->SetPoint(1,0.9198646,0.228013);
699       graph->SetPoint(2,0.996614,0.5895765);
700       graph->SetPoint(3,1.041761,0.8241042);
701       graph->SetPoint(4,1.059819,0.8794788);
702       graph->SetPoint(5,1.077878,0.9348534);
703       graph->SetPoint(6,1.100451,0.980456);
704       graph->SetPoint(7,1.141084,0.9967427);
705       graph->SetPoint(8,1.204289,0.9153094);
706       graph->SetPoint(9,1.276524,0.742671);
707       graph->SetPoint(10,1.402935,0.465798);
708       graph->SetPoint(11,1.515801,0.3029316);
709       graph->SetPoint(12,1.73702,0.1465798);
710       graph->SetPoint(13,1.985327,0.08143322);
711       graph->SetPoint(14,2.301354,0.04234528);
712       graph->SetPoint(15,2.56772,0.02931596);
713 #endif
714     }
715     return graph;
716   }
717
718   /** Get the correction to Bethe-Bloc from Review of Particle Physics
719       (fig 27.8). 
720   */
721   TGraph* GetCorr() 
722   {
723     static TGraph* graph = 0;
724     if (!graph) {
725       graph = new TGraph(14);
726       graph->SetName("graph");
727       graph->SetTitle("(#Delta_{p}/x)/(dE/dx)|_{mip} for 320#mu Si");
728       graph->GetHistogram()->SetXTitle("#beta#gamma = p/m");
729       graph->SetFillColor(1);
730       graph->SetLineColor(7);
731       graph->SetMarkerStyle(21);
732       graph->SetMarkerSize(0.6);
733       graph->SetPoint(0,1.196058,0.9944915);
734       graph->SetPoint(1,1.28502,0.9411017);
735       graph->SetPoint(2,1.484334,0.8559322);
736       graph->SetPoint(3,1.984617,0.7491525);
737       graph->SetPoint(4,2.658367,0.6983051);
738       graph->SetPoint(5,3.780227,0.6779661);
739       graph->SetPoint(6,4.997358,0.6741525);
740       graph->SetPoint(7,8.611026,0.684322);
741       graph->SetPoint(8,15.28296,0.6995763);
742       graph->SetPoint(9,41.54516,0.7186441);
743       graph->SetPoint(10,98.91461,0.7288136);
744       graph->SetPoint(11,203.2734,0.7326271);
745       graph->SetPoint(12,505.6421,0.7338983);
746       graph->SetPoint(13,896.973,0.7338983);
747     }
748     return graph;
749   }
750   
751   ClassDef(DrawHits,0);
752 };
753
754 //____________________________________________________________________
755 //
756 // EOF
757 //