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