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 //