AliRoot Core  3dc7879 (3dc7879)
DrawHits.C
Go to the documentation of this file.
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 
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  {
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 
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  }
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 
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 
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 
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 
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 
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 
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 //
DrawHits(const char *pdgName="pi+", Int_t m=500, Double_t emin=1, Double_t emax=1000, Int_t n=900, Double_t tmin=1e-2, Double_t tmax=1e3)
Definition: DrawHits.C:52
Draw hit energy loss.
Definition: DrawHits.C:41
Float_t M() const
Definition: AliFMDHit.cxx:165
TGraph * FromRPPMean()
Definition: DrawHits.C:477
TStyle * gStyle
Hit in the FMD.
FMD utility classes for reading FMD data.
TGraph * FromRPPFull()
Definition: DrawHits.C:354
Float_t P() const
Definition: AliFMDHit.cxx:156
TH2D * fElossVsPMQ
Definition: DrawHits.C:44
Float_t p[]
Definition: kNNTest.C:133
TGraph * FromRPPNoDelta()
Definition: DrawHits.C:396
Int_t Track() const
Definition: AliHit.h:24
TH1D * fBetaGamma
Definition: DrawHits.C:46
void ScaleGraph(TGraph *graph, bool density=true, double mass=1)
Definition: DrawHits.C:344
Bool_t IsPhysicalPrimary(Int_t i, Bool_t useInEmbedding=kFALSE)
Definition: AliStack.cxx:1047
Double_t chi2
Definition: AnalyzeLaser.C:7
Float_t Edep() const
Definition: AliFMDHit.h:82
TGraph * GetResp()
Definition: DrawHits.C:631
const char * GetName() const
Definition: AliFMDHit.cxx:133
Float_t Q() const
Definition: AliFMDHit.cxx:175
TGraph * FromGFMATE()
Definition: DrawHits.C:514
static TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max)
void ShowFit(Double_t x1, Double_t y1, const char *title, TF1 *f, Double_t dx=0, Double_t dy=0.05)
Definition: DrawHits.C:148
virtual void AddLoad(ETrees tree)
Definition: AliFMDInput.h:134
TParticlePDG * fPdg
Definition: DrawHits.C:47
AliStack * fStack
Definition: AliFMDInput.h:416
TLatex latex
TF1 * f
Definition: interpolTest.C:21
TH1D * fEloss
Definition: DrawHits.C:45
Base class for reading in various FMD data. The class loops over all found events. For each event the specified data is read in. The class then loops over all elements of the read data, and process these with user defined code.
Definition: AliFMDInput.h:106
const Double_t fBetaGammaMip
Definition: DrawHits.C:49
Float_t Length() const
Definition: AliFMDHit.h:100
Bool_t Finish()
Definition: DrawHits.C:178
AliFMDhit is the hit class for the FMD. Hits are the information that comes from a Monte Carlo at eac...
Definition: AliFMDHit.h:30
Int_t fEventCount
Definition: AliFMDInput.h:444
void res(Char_t i)
Definition: Resolution.C:2
Bool_t IsStop() const
Definition: AliFMDHit.h:102
const Double_t fRho
Definition: DrawHits.C:48
Bool_t ProcessHit(AliFMDHit *hit, TParticle *p)
Definition: DrawHits.C:95
TGraph * GetCorr()
Definition: DrawHits.C:721
TGraph * FromRPPNoRad()
Definition: DrawHits.C:438