FMD/scripts/DrawESD.C

   1 //____________________________________________________________________
   2 //
   3 // $Id$
   4 //
   5 // Script that contains a class to draw eloss from hits, versus ADC
   6 // counts from digits, using the AliFMDInputHits class in the util library.
   7 //
   8 // It draws the energy loss versus the p/(mq^2).  It can be overlayed
   9 // with the Bethe-Bloc curve to show how the simulation behaves
  10 // relative to the expected.
  11 //
  12 // Use the script `Compile.C' to compile this class using ACLic.
  13 //
  14 #include <TH1D.h>
  15 #include <AliFMDHit.h>
  16 #include <AliFMDDigit.h>
  17 #include <AliFMDInput.h>
  18 #include <AliFMDUShortMap.h>
  19 #include <AliFMDFloatMap.h>
  20 #include <AliFMDRecPoint.h>
  21 #include <AliESDFMD.h>
  22 #include <AliLog.h>
  23 #include <iostream>
  24 #include <TStyle.h>
  25 #include <TArrayF.h>
  26 #include <TCanvas.h>
  27 #include <TMath.h>
  28 #include <TF1.h>
  29 #include <TSpectrum.h>
  30 #include <TLegend.h>
  31 #include <TLatex.h>
  32 #include <TLine.h>
  33 #include <TPolyMarker.h>
  34 #include <TSpectrum.h>
  35 #include <TList.h>
  36 #include <TGraph.h>
  37
  38 Double_t landau(Double_t* xp, Double_t* pp)
  39 {
  40   Double_t x     = xp[0];
  41   Double_t A     = pp[0];
  42   Double_t mpv   = pp[1];
  43   Double_t w     = pp[2];
  44
  45   Double_t v     = mpv; //  + w * 0.22278;
  46   return A * TMath::Landau(x, v, w);
  47 }
  48
  49 Double_t foldLandau(Double_t* xp, Double_t* pp)
  50 {
  51   Double_t x     = xp[0];
  52   Double_t A     = pp[0];
  53   Double_t mpv   = pp[1];
  54   Double_t w     = pp[2];
  55   Double_t B     = pp[3];
  56   Double_t sigma = pp[4];
  57
  58   return A * (TMath::Landau(x,mpv,w) + B * TMath::Gaus(x, mpv, sigma));
  59 }
  60
  61 /** @class DrawESD
  62     @brief Draw digit ADC versus Rec point mult
  63     @code
  64     Root> .L Compile.C
  65     Root> Compile("DrawESD.C")
  66     Root> DrawESD c
  67     Root> c.Run();
  68     @endcode
  69     @ingroup FMD_script
  70  */
  71 class DrawESD : public AliFMDInput
  72 {
  73 private:
  74   TH1D* fMult; // Histogram
  75   const Double_t fCorr;
  76   TList fCleanup;
  77 public:
  78   //__________________________________________________________________
  79   DrawESD(Int_t n=2000, Double_t mmin=-0.5, Double_t mmax=15.5)
  80     : fCorr(1) // 0.68377 / 1.1)
  81   {
  82     AddLoad(kESD);
  83     fMult = new TH1D("mult", "#DeltaE/#DeltaE_{MIP)", n, mmin, mmax);
  84     fMult->Sumw2();
  85     fMult->SetFillColor(kRed+1);
  86     fMult->SetFillStyle(3001);
  87     fMult->SetXTitle("#DeltaE/#DeltaE_{MIP}");
  88     fCleanup.Add(fMult);
  89   }
  90   ~DrawESD()
  91   {
  92     fCleanup.Delete();
  93   }
  94   //__________________________________________________________________
  95   /** Begining of event
  96       @param ev Event number
  97       @return @c false on error */
  98   Bool_t Begin(Int_t ev)
  99   {
 100     return AliFMDInput::Begin(ev);
 101   }
 102   //__________________________________________________________________
 103   Bool_t ProcessESD(UShort_t /* det */,
 104                     Char_t   /* rng */,
 105                     UShort_t /* sec */,
 106                     UShort_t /* str */,
 107                     Float_t  eta,
 108                     Float_t  mult)
 109   {
 110     // Cache the energy loss
 111     // if (mult > 0)
 112     //   Info("ProcessESD", "FMD%d%c[%2d,%3d]=%f", det, rng, sec, str, mult);
 113     if (mult <= 0 || mult == AliESDFMD::kInvalidMult) return kTRUE;
 114     Double_t x = mult;
 115     if (!fESD->IsAngleCorrected()) {
 116       Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
 117       Double_t corr  = TMath::Abs(TMath::Cos(theta));
 118       Double_t cmult = corr * mult;
 119       x = cmult;
 120     }
 121     if (x > 0.001) fMult->Fill(x);
 122     return kTRUE;
 123   }
 124   //__________________________________________________________________
 125   Bool_t ProcessESDs()
 126   {
 127     if (!AliFMDInput::ProcessESDs()) return kFALSE;
 128     // Info("ProcessESDs", "ESD is %sangle corrected",
 129     //      fESD->IsAngleCorrected() ? "" : "not ");
 130     return kTRUE;
 131   }
 132   //__________________________________________________________________
 133   TF1* FitPeak(Int_t n, TH1D* hist, Double_t min, Double_t& max)
 134   {
 135     if (TMath::Abs(max-min) < .25) return 0;
 136     std::cout << "Fit peak in range " << min << " to " << max << std::endl;
 137     TF1* l = new TF1(Form("l%02d", n), "landau", min, max);
 138     hist->GetXaxis()->SetRangeUser(0, 4);
 139     hist->Fit(l, "0Q", "", min, max);
 140     Double_t mpv   = l->GetParameter(1);
 141     Double_t empv  = l->GetParError(1);
 142     Double_t sigma = l->GetParameter(2);
 143     l->SetRange(mpv-empv, mpv+3*sigma);
 144     hist->Fit(l, "EMQ0", "", mpv-3*empv, mpv+3*sigma);
 145     std::cout << "Peak # " << n << " [" << min << "," << max << "]\n"
 146               << " MPV: " << l->GetParameter(1)
 147               << " +/- "  << l->GetParError(1)
 148               << " Var: " << l->GetParameter(2)
 149               << " +/- "  << l->GetParError(2)
 150               << " Chi^2/NDF: " << l->GetChisquare() / l->GetNDF()
 151               << std::endl;
 152     mpv   = l->GetParameter(1);
 153     sigma = l->GetParameter(2);
 154     min   = mpv - sigma * 2; // (n==1 ? 3 : 2);
 155     max   = mpv + sigma * 3;
 156     // l->SetRange(min, max);
 157     l->Draw("same");
 158     return  l;
 159   }
 160   //__________________________________________________________________
 161   void MaxInRange(TH1D* hist, Double_t min, Double_t& mean, Double_t& var)
 162   {
 163     hist->GetXaxis()->SetRangeUser(min, 4);
 164     mean = hist->GetMean();
 165     var  = hist->GetRMS();
 166   }
 167
 168   //__________________________________________________________________
 169   const char* PrettyFloat(float x)
 170   {
 171     if (x == 0) return Form("%9.4f", x);
 172     float e = TMath::Floor(TMath::Log10(TMath::Abs(x)));
 173     if (TMath::Abs(e) < 4) {
 174       return Form("%9.4f", x);
 175     }
 176     float f = x / TMath::Power(10,e);
 177     return Form("%4.2f#times10^{%d}", f, int(e));
 178   }
 179   //__________________________________________________________________
 180   void ShowFit(Double_t x1, Double_t y1, const char* title,
 181                TF1* f, Double_t dx=0, Double_t dy=0.05)
 182   {
 183     Double_t x = x1, y = y1;
 184     TLatex* latex = new TLatex(x, y, title);
 185     latex->SetTextFont(132);
 186     latex->SetTextSize(0.8*dy);
 187     latex->SetNDC();
 188     latex->Draw();
 189     x -= dx;
 190     y -= dy;
 191     const Double_t eqDx=0.1;
 192     Double_t chi2 = f->GetChisquare();
 193     Int_t    ndf  = f->GetNDF();
 194     Double_t prob = f->GetProb();
 195     latex->DrawLatex(x, y, "#chi^{2}/NDF");
 196     latex->DrawLatex(x+eqDx, y, Form("= %7.4f/%3d=%5.2f (%7.3f%%)",
 197                                      chi2, ndf, chi2/ndf, 100*prob));
 198     Int_t     n = f->GetNpar();
 199     Double_t* p = f->GetParameters();
 200     Double_t* e = f->GetParErrors();
 201     for (int i = 0; i < n; i++) {
 202       x -= dx;
 203       y -= dy;
 204       latex->DrawLatex(x, y, f->GetParName(i));
 205       latex->DrawLatex(x+eqDx, y, Form("= %s", PrettyFloat(p[i])));
 206       latex->DrawLatex(x+2.2*eqDx, y, Form("#pm %s", PrettyFloat(e[i])));
 207     }
 208   }
 209   //__________________________________________________________________
 210   Bool_t Finish()
 211   {
 212     Info("Finish", "Will draw results");
 213     gStyle->SetPalette(1);
 214     gStyle->SetOptTitle(0);
 215     gStyle->SetOptFit(1111111);
 216     gStyle->SetCanvasColor(0);
 217     gStyle->SetCanvasBorderSize(0);
 218     gStyle->SetPadColor(0);
 219     gStyle->SetPadBorderSize(0);
 220
 221     if (fMult->GetEntries() <= 0) return kFALSE;
 222
 223     TCanvas* c = new TCanvas("c", "C", 1200, 800);
 224     fCleanup.Add(c);
 225     c->cd();
 226     c->SetLogy();
 227     c->SetTopMargin(0.05);
 228     c->SetRightMargin(0.05);
 229     c->SetFillColor(0);
 230     c->SetBorderMode(0);
 231
 232     TLegend* leg = new TLegend(.1, .1, .4, .2);
 233     leg->SetFillColor(0);
 234     leg->SetBorderSize(1);
 235
 236     DrawMult(c, leg);
 237
 238     Double_t xmax=0, xmin=0, ymax=0;
 239     FindMinMax(xmin, xmax, ymax);
 240
 241     FitLandau(xmin, xmax, ymax, leg);
 242     DrawResponse(xmax, ymax, leg);
 243
 244     // TF1* f = FitMultiLandau(xmin, xmax, leg);
 245     // DrawLandaus(f);
 246
 247     c->cd();
 248     leg->Draw();
 249
 250     c->Modified();
 251     c->Update();
 252     c->cd();
 253     c->SaveAs("esd_eloss.png");
 254
 255     fCleanup.Add(leg);
 256
 257     return kTRUE;
 258
 259   }
 260   //__________________________________________________________________
 261   /**
 262    * Draw the multiplicity distribution
 263    *
 264    */
 265   void DrawMult(TCanvas* /* c */, TLegend* leg)
 266   {
 267     fMult->GetXaxis()->SetRangeUser(0.2,20);
 268     Double_t integral = fMult->Integral();
 269     Info("DrawMult", "Integral in range [0.2,20]=%f (%f entries)",
 270          integral, fMult->GetEntries());
 271     fMult->Scale(1. / integral);
 272     Info("DrawMult", "Integral in range [0.2,20]=%f (%f entries)",
 273          fMult->Integral(), fMult->GetEntries());
 274     Double_t max = 1.5 * fMult->GetMaximum();
 275     fMult->GetXaxis()->SetRangeUser(0,4);
 276     fMult->SetMaximum(max);
 277     fMult->SetStats(kFALSE);
 278     fMult->Draw("e same");
 279     leg->AddEntry(fMult, "Strip signal", "lf");
 280
 281   }
 282   //__________________________________________________________________
 283   /**
 284    * Find the minimum and maximum values in range
 285    *
 286    * @param xmin
 287    * @param xmax
 288    * @param ymax
 289    */
 290   void FindMinMax(Double_t& xmin, Double_t& xmax, Double_t& ymax)
 291   {
 292     fMult->GetXaxis()->SetRangeUser(0.1,4);
 293     TSpectrum    s(10);
 294     Int_t        nPeak  = s.Search(fMult);
 295     fMult->GetXaxis()->SetRangeUser(0.1, 4);
 296     Int_t        bmax   = fMult->GetMaximumBin();
 297     xmax                = fMult->GetBinCenter(bmax);
 298     fMult->GetXaxis()->SetRangeUser(0.1, xmax);
 299     Double_t     bmin   = fMult->GetMinimumBin();
 300     xmin                = fMult->GetBinCenter(bmin);
 301     ymax                = fMult->GetBinContent(bmax);
 302     Info("Finish", "Simple peak finding found x_max=%f, x_min=%f y_max=%g",
 303          xmax, xmin, ymax);
 304
 305     if (nPeak > 1) {
 306       TPolyMarker* pm = static_cast<TPolyMarker*>(fMult->GetListOfFunctions()
 307                                                   ->FindObject("TPolyMarker"));
 308
 309       // Peaks are ordered by size
 310       Double_t* peakX = pm->GetX();
 311       Double_t* peakY = pm->GetY();
 312       Double_t  xlmax = peakX[1];
 313       xmax            = peakX[0];
 314       ymax            = peakY[0];
 315       if (xmax < xlmax) {
 316         xmax  = xlmax;
 317         xlmax = peakX[0];
 318         ymax  = peakY[1];
 319       }
 320
 321       fMult->GetXaxis()->SetRangeUser(xlmax, xmax);
 322       bmin  = fMult->GetMinimumBin();
 323       xmin  = fMult->GetBinCenter(bmin);
 324       Info("Finish", "Spectrum peak finding found x_max=%f, x_min=%f y_max=%g",
 325            xmax, xmin, ymax);
 326     }
 327   }
 328   //__________________________________________________________________
 329   /**
 330    * Fit a landau and a landau+gaussian to the data
 331    *
 332    * @param xmin  Minimum of peak range
 333    * @param xmax  Maximum of peak range
 334    * @param ymax  Y value in MIP peak
 335    * @param leg   Legend
 336    */
 337   void FitLandau(Double_t xmin, Double_t xmax, Double_t& ymax, TLegend* leg)
 338   {
 339     fMult->GetXaxis()->SetRangeUser(xmin, xmax+(xmax-xmin));
 340     Double_t mean = fMult->GetMean();
 341     Double_t var  = fMult->GetRMS();
 342     Info("Finish", "Peak range [%f,%f] mean=%f, var=%f",
 343          xmin, xmax+(xmax-xmin), mean, var);
 344     Double_t lowCut = mean-var;
 345     Double_t hiCut  = 4; // 2*mean;
 346     Info("Finish", "Low cut set to %f", lowCut);
 347     fMult->GetXaxis()->SetRangeUser(0, hiCut);
 348
 349     TF1* pl = MakeLandau(lowCut, hiCut, xmax, var/2);
 350     fMult->Fit(pl, "NEM", "", lowCut, hiCut);
 351     ymax = pl->GetMaximum();
 352     Info("Finish", "Maximum of landau is at %f (y_max=%g)",
 353          pl->GetMaximumX(), ymax);
 354
 355     TF1* gl  = MakeFoldLandau(lowCut, hiCut, pl, xmax, var/2);
 356     gl->SetLineColor(kRed+1);
 357     // gl->SetParLimits(1, xmax-var, xmax+var);
 358     fMult->Fit(gl, "NEM", "", lowCut, hiCut);
 359     TF1* l  = MakeLandau(lowCut, hiCut, xmax, var/2);
 360     l->SetLineColor(kGreen+1);
 361     l->SetParameters(gl->GetParameter(0),
 362                      gl->GetParameter(1),
 363                      gl->GetParameter(2));
 364     TF1* g  = new TF1("g", "gaus", lowCut, hiCut);
 365     g->SetParNames("A", "#mu", "#sigma");
 366     g->SetLineColor(kBlue+1);
 367     g->SetParameters(gl->GetParameter(3)*gl->GetParameter(0),
 368                      gl->GetParameter(1),
 369                      gl->GetParameter(2));
 370     fMult->GetXaxis()->SetRangeUser(0,4);
 371     fMult->DrawCopy("E");
 372     fMult->DrawCopy("HIST SAME");
 373     pl->Draw("same");
 374     gl->Draw("same");
 375     g->Draw("Same");
 376     l->Draw("Same");
 377     fCleanup.Add(pl);
 378     fCleanup.Add(l);
 379     fCleanup.Add(g);
 380     fCleanup.Add(gl);
 381
 382     ShowFit(.5, .90, "Landau", pl);
 383     ShowFit(.5, .65, "Landau+Gaussian", gl);
 384
 385     leg->AddEntry(pl, Form("Landau fit - #chi^{2}/NDF=%f",
 386                            pl->GetChisquare()/pl->GetNDF()), "l");
 387     leg->AddEntry(gl, Form("Landau+Gaussian fit - #chi^{2}/NDF=%f",
 388                            gl->GetChisquare()/gl->GetNDF()), "l");
 389     leg->AddEntry(l, "Landau part", "l");
 390     leg->AddEntry(g, "Gaussian part", "l");
 391   }
 392
 393   //__________________________________________________________________
 394   /**
 395    * Superimpose the response graph from the RPP
 396    *
 397    * @param ymax Y value of multiplicity spectra in the landau peak
 398    * @param leg  Legend
 399    */
 400   void DrawResponse(Double_t xmax, Double_t ymax, TLegend* leg)
 401   {
 402     TGraph*   resp = GetResp();
 403     // TGraph*   corr = GetCorr();
 404     Double_t* x    = resp->GetX();
 405     Double_t* y    = resp->GetY(); // [MeV cm^2/g]
 406     TGraph*   gr   = new TGraph(resp->GetN());
 407     gr->SetName(Form("%sCorr", resp->GetName()));
 408     gr->SetTitle(resp->GetTitle());
 409     gr->SetLineStyle(resp->GetLineStyle());
 410     gr->SetLineColor(kMagenta+1);
 411     gr->SetLineWidth(2);
 412     TGraph*   gr2   = new TGraph(resp->GetN());
 413     gr2->SetName(Form("%sCorr", resp->GetName()));
 414     gr2->SetTitle(resp->GetTitle());
 415     gr2->SetLineStyle(resp->GetLineStyle());
 416     gr2->SetLineColor(kCyan+1);
 417     gr2->SetLineWidth(2);
 418     // const Double_t rho = 2.33;   // [g/cm^3]
 419     // const Double_t thk = 0.320;  // [cm]
 420     const Double_t mip = 1.664;  // [MeV cm^2/g]
 421     // const Double_t bgm = 3.4601; // beta*gamma of a MIP
 422     // Double_t  xs2 = corr->Eval(bgm); // [1]
 423     // Double_t  xss = 1.1;
 424     Double_t  xs  = 1/mip;
 425     Double_t xxmax = 0;
 426     Double_t yymax = 0;
 427     for (Int_t i = 0; i < gr->GetN(); i++) {
 428       if (y[i] > yymax) {
 429         yymax = y[i];
 430         xxmax = xs * x[i];
 431       }
 432       gr->SetPoint(i, x[i] * xs, y[i] * ymax);
 433     }
 434     Info("DrawResponse", "Maximum at x=%f (xmax=%f)", xxmax, xmax);
 435     Double_t xs2 = xmax / xxmax;
 436     Info("DrawResponse", "Correction factor: %f", xs2);
 437     for (Int_t i = 0; i < gr->GetN(); i++) {
 438       gr2->SetPoint(i, x[i] * xs * xs2, y[i] * ymax);
 439     }
 440     gr->Draw("C same");
 441     gr2->Draw("C same");
 442
 443     leg->AddEntry(gr, "Response", "l");
 444     leg->AddEntry(gr2, "Response", "l");
 445   }
 446   //__________________________________________________________________
 447   /**
 448    * Fit sum of landaus to the multiplicity distribution
 449    *
 450    * @param xmin Minimum of range
 451    * @param xmax Maximum of range
 452    * @param leg  Legend
 453    *
 454    * @return Fitted function
 455    */
 456   TF1* FitMultiLandau(Double_t xmin, Double_t xmax, TLegend* leg)
 457   {
 458     fMult->GetXaxis()->SetRangeUser(xmin, xmax+(xmax-xmin));
 459     Double_t mean = fMult->GetMean();
 460     Double_t rms  = fMult->GetRMS();
 461
 462     Double_t x1   = mean-rms/2; // .75; // .8;  // .65 / fCorr;
 463     Double_t x2   = mean+rms/2; // 1.3; // 1.7; // fCorr;
 464     TF1*     l1   = FitPeak(1, fMult, x1, x2);
 465     x2            = TMath::Max(mean+rms/2, x2);
 466     MaxInRange(fMult, x2, mean, rms);
 467     Double_t x3   = mean + rms;
 468     TF1*     l2   = FitPeak(2, fMult, x2, x3);
 469     TF1*     f    = 0;
 470     Double_t diff = 0;
 471     if (l2) {
 472       diff          = l2->GetParameter(1)-l1->GetParameter(1);
 473       f             = new TF1("user", "landau(0)+landau(3)", x1, x3);
 474       f->SetParNames("A_{1}", "Mpv_{1}", "#sigma_{1}",
 475                      "A_{2}", "Mpv_{2}", "#sigma_{2}",
 476                      "A_{3}", "Mpv_{3}", "#sigma_{3}");
 477       f->SetParameters(l1->GetParameter(0),
 478                        l1->GetParameter(1),
 479                        l1->GetParameter(2),
 480                        l2 ? l2->GetParameter(0) : 0,
 481                        l2 ? l2->GetParameter(1) : 0,
 482                        l2 ? l2->GetParameter(2) : 0,
 483                        l2->GetParameter(0)/10,
 484                        l2->GetParameter(1) + diff,
 485                        l2->GetParameter(2));
 486     }
 487     else {
 488       x3 = x2;
 489       f  = new TF1("usr", "landau", x1, x3);
 490     }
 491
 492     std::cout << "Range: " << x1 << "-" << x3 << std::endl;
 493
 494     fMult->GetXaxis()->SetRangeUser(0, 4);
 495     fMult->Fit(f, "NQR", "", x1, x3);
 496     fMult->Fit(f, "NMER", "E1", x1, x3);
 497
 498     l1->SetLineColor(3);
 499     l1->SetLineWidth(2);
 500     l1->SetRange(0, 4);
 501     l1->Draw("same");
 502     if (l2) {
 503       l2->SetLineColor(4);
 504       l2->SetLineWidth(2);
 505       l2->SetRange(0, 4);
 506       l2->Draw("same");
 507     }
 508     f->SetLineWidth(2);
 509     f->SetRange(0, 4);
 510     f->Draw("same");
 511
 512     fCleanup.Add(l1);
 513     if (l2) fCleanup.Add(l2);
 514     fCleanup.Add(f);
 515
 516
 517     leg->AddEntry(l1, "1 particle Landau", "l");
 518     if (l2) leg->AddEntry(l2, "2 particle Landau", "l");
 519     leg->AddEntry(f,  "1+2 particle Landau", "l");
 520
 521     return f;
 522   }
 523   //__________________________________________________________________
 524   /**
 525    * Draw landau functions in a separate canvas
 526    *
 527    * @param f Multi-landau function
 528    */
 529   void DrawLandaus(TF1* f)
 530   {
 531     TCanvas* c = new TCanvas("c2", "Landaus");
 532     // c->SetLogy();
 533     fMult->DrawClone("axis");
 534     f->Draw("same");
 535     Double_t* p1 = f->GetParameters();
 536     Double_t* p2 = &(p1[3]);
 537     TF1* ll1 = new TF1("ll1", "landau", 0, 4);
 538     ll1->SetParameters(p1);
 539     ll1->SetLineColor(3);
 540     ll1->Draw("same");
 541     TF1* ll2 = new TF1("ll2", "landau", 0, 4);
 542     ll2->SetParameters(p2);
 543     ll2->SetLineColor(4);
 544     ll2->Draw("same");
 545
 546     Double_t diff = ll2->GetParameter(1)-ll1->GetParameter(1);
 547     Double_t y1   = fMult->GetMinimum() * 1.1;
 548     Double_t y2   = fMult->GetMaximum() * .9;
 549     Double_t xc1  = p1[1]-3*p1[2];
 550     Double_t xc2  = p2[1]-2*p2[2];
 551     Double_t xc3  = p2[1]-2*p2[2]+diff;
 552     TLine* c1 = new TLine(xc1, y1, xc1, y2);
 553     c1->Draw("same");
 554     TLine* c2 = new TLine(xc2, y1, xc2, y2);
 555     c2->Draw("same");
 556     TLine* c3 = new TLine(xc3, y1, xc3, y2);
 557     c3->Draw("same");
 558
 559     TLegend* l = new TLegend(0.6, 0.6, .89, .89);
 560     l->AddEntry(ll1, "1 particle Landau", "l");
 561     l->AddEntry(ll2, "2 particle Landau", "l");
 562     l->AddEntry(f,  "1+2 particle Landau", "l");
 563     l->SetFillColor(0);
 564     l->Draw("same");
 565
 566     c->Modified();
 567     c->Update();
 568     c->cd();
 569   }
 570
 571   //__________________________________________________________________
 572   /**
 573    * Get the response functin @f$ f(\Delta_p/x)@f$ from Review of
 574    * Particle Physics (fig. 27.7).  It is scaled to the value at MPV.
 575    *
 576    * @return Graph of response
 577    */
 578   TGraph* GetResp()
 579   {
 580     static TGraph*  graph = 0;
 581     if (!graph) {
 582       graph = new TGraph;
 583       graph->SetName("si_resp");
 584       graph->SetTitle("f(#Delta/x) scaled to the MPV value ");
 585       graph->GetHistogram()->SetXTitle("#Delta/x (MeVcm^{2}/g)");
 586       graph->GetHistogram()->SetYTitle("f(#Delta/x)");
 587       graph->SetLineColor(kBlue+1);
 588       graph->SetLineWidth(2);
 589       graph->SetFillColor(kBlue+1);
 590       graph->SetMarkerStyle(21);
 591       graph->SetMarkerSize(0.6);
 592 #if 0
 593       // Figure 27.7 or Review of Particle physics - Straggeling function in
 594       // silicon of 500 MeV Pions, normalised to unity at the most probable
 595       // value.
 596       graph->SetPoint(0,0.808094,0.00377358);
 597       graph->SetPoint(1,0.860313,0.0566038);
 598       graph->SetPoint(2,0.891645,0.116981);
 599       graph->SetPoint(3,0.912533,0.181132);
 600       graph->SetPoint(4,0.928198,0.260377);
 601       graph->SetPoint(5,0.938642,0.320755);
 602       graph->SetPoint(6,0.954308,0.377358);
 603       graph->SetPoint(7,0.964752,0.433962);
 604       graph->SetPoint(8,0.975196,0.490566);
 605       graph->SetPoint(9,0.98564,0.550943);
 606       graph->SetPoint(10,0.996084,0.611321);
 607       graph->SetPoint(11,1.00653,0.667925);
 608       graph->SetPoint(12,1.02219,0.732075);
 609       graph->SetPoint(13,1.03264,0.784906);
 610       graph->SetPoint(14,1.0483,0.845283);
 611       graph->SetPoint(15,1.06397,0.901887);
 612       graph->SetPoint(16,1.09008,0.958491);
 613       graph->SetPoint(17,1.10574,0.984906);
 614       graph->SetPoint(18,1.13708,1);
 615       graph->SetPoint(19,1.13708,1);
 616       graph->SetPoint(20,1.15796,0.988679);
 617       graph->SetPoint(21,1.17363,0.966038);
 618       graph->SetPoint(22,1.19974,0.916981);
 619       graph->SetPoint(23,1.2154,0.89434);
 620       graph->SetPoint(24,1.23629,0.837736);
 621       graph->SetPoint(25,1.2624,0.784906);
 622       graph->SetPoint(26,1.28329,0.724528);
 623       graph->SetPoint(27,1.3094,0.664151);
 624       graph->SetPoint(28,1.32507,0.611321);
 625       graph->SetPoint(29,1.3564,0.550943);
 626       graph->SetPoint(30,1.41384,0.445283);
 627       graph->SetPoint(31,1.44517,0.392453);
 628       graph->SetPoint(32,1.48695,0.335849);
 629       graph->SetPoint(33,1.52872,0.286792);
 630       graph->SetPoint(34,1.58094,0.237736);
 631       graph->SetPoint(35,1.63838,0.196226);
 632       graph->SetPoint(36,1.68016,0.169811);
 633       graph->SetPoint(37,1.75326,0.135849);
 634       graph->SetPoint(38,1.81593,0.113208);
 635       graph->SetPoint(39,1.89426,0.0981132);
 636       graph->SetPoint(40,1.96214,0.0830189);
 637       graph->SetPoint(41,2.0718,0.0641509);
 638       graph->SetPoint(42,2.19191,0.0490566);
 639       graph->SetPoint(43,2.31723,0.0415094);
 640       graph->SetPoint(44,2.453,0.0301887);
 641       graph->SetPoint(45,2.53133,0.0264151);
 642       graph->SetPoint(46,2.57833,0.0264151);
 643 #else
 644       graph->SetPoint(0,0.8115124,0.009771987);
 645       graph->SetPoint(1,0.9198646,0.228013);
 646       graph->SetPoint(2,0.996614,0.5895765);
 647       graph->SetPoint(3,1.041761,0.8241042);
 648       graph->SetPoint(4,1.059819,0.8794788);
 649       graph->SetPoint(5,1.077878,0.9348534);
 650       graph->SetPoint(6,1.100451,0.980456);
 651       graph->SetPoint(7,1.141084,0.9967427);
 652       graph->SetPoint(8,1.204289,0.9153094);
 653       graph->SetPoint(9,1.276524,0.742671);
 654       graph->SetPoint(10,1.402935,0.465798);
 655       graph->SetPoint(11,1.515801,0.3029316);
 656       graph->SetPoint(12,1.73702,0.1465798);
 657       graph->SetPoint(13,1.985327,0.08143322);
 658       graph->SetPoint(14,2.301354,0.04234528);
 659       graph->SetPoint(15,2.56772,0.02931596);
 660 #endif
 661     }
 662     return graph;
 663   }
 664   //__________________________________________________________________
 665   /**
 666    * Get the correction to Bethe-Bloc from Review of Particle Physics
 667    * (fig 27.8).
 668    *
 669    * @return correction graph
 670    */
 671   TGraph* GetCorr()
 672   {
 673     static TGraph* graph = 0;
 674     if (!graph) {
 675       graph = new TGraph(14);
 676       graph->SetName("graph");
 677       graph->SetTitle("(#Delta_{p}/x)/(dE/dx)|_{mip} for 320#mu Si");
 678       graph->GetHistogram()->SetXTitle("#beta#gamma = p/m");
 679       graph->GetHistogram()->SetYTitle("#frac{#Delta_{p}/x)}{(dE/dx)|_{mip}}");
 680       graph->SetFillColor(1);
 681       graph->SetLineColor(7);
 682       graph->SetMarkerStyle(21);
 683       graph->SetMarkerSize(0.6);
 684       graph->SetPoint(0,1.196058,0.9944915);
 685       graph->SetPoint(1,1.28502,0.9411017);
 686       graph->SetPoint(2,1.484334,0.8559322);
 687       graph->SetPoint(3,1.984617,0.7491525);
 688       graph->SetPoint(4,2.658367,0.6983051);
 689       graph->SetPoint(5,3.780227,0.6779661);
 690       graph->SetPoint(6,4.997358,0.6741525);
 691       graph->SetPoint(7,8.611026,0.684322);
 692       graph->SetPoint(8,15.28296,0.6995763);
 693       graph->SetPoint(9,41.54516,0.7186441);
 694       graph->SetPoint(10,98.91461,0.7288136);
 695       graph->SetPoint(11,203.2734,0.7326271);
 696       graph->SetPoint(12,505.6421,0.7338983);
 697       graph->SetPoint(13,896.973,0.7338983);
 698     }
 699     return graph;
 700   }
 701   //__________________________________________________________________
 702   /**
 703    * Make a Landau function object.
 704    *
 705    * @param min  Minimum of fit range
 706    * @param max  Maximum of fit range
 707    * @param p    Peak position
 708    * @param v    Variance around peak
 709    *
 710    * @return Landau function object
 711    */
 712   TF1* MakeLandau(Double_t min, Double_t max, Double_t p, Double_t v)
 713   {
 714     TF1* f = new TF1("l", "landau", min, max);
 715     f->SetParNames("A", "#delta", "#xi");
 716     f->SetParameters(1, p, p/10);
 717     if      (false) f->FixParameter(1,p);
 718     else if (false) f->SetParLimits(1, p-v, p+v);
 719     return f;
 720   }
 721   //__________________________________________________________________
 722   /**
 723    * Make a Landau, folded with a gaussian, function object
 724    *
 725    * @param min  Minimum of fit range
 726    * @param max  Maximum of fit range
 727    * @param l    Seed Landau function object
 728    * @param p    Peak position
 729    * @param v    Variance around peak
 730    *
 731    * @return Landau+Gaus function object
 732    */
 733   TF1* MakeFoldLandau(Double_t min, Double_t max, TF1* l,
 734                       Double_t p, Double_t v)
 735   {
 736     TF1* f = new TF1("gl", &foldLandau, min, max, 5);
 737     f->SetParNames("A", "#delta", "#xi", "B", "#sigma");
 738     f->SetParameters(l->GetParameter(0),
 739                      l->GetParameter(1),
 740                      l->GetParameter(2),
 741                      l->GetParameter(0)/1000,
 742                      l->GetParameter(2));
 743     f->SetParLimits(3, 1e-7, 1e-1);
 744     if      (false) f->FixParameter(1,p);
 745     else if (true)  f->SetParLimits(1, p-2*v, p+2*v);
 746     return f;
 747   }
 748
 749   ClassDef(DrawESD,0);
 750
 751 };
 752
 753 //____________________________________________________________________
 754 //
 755 // EOF
 756 //