]>
Commit | Line | Data |
---|---|---|
bf000c32 | 1 | //____________________________________________________________________ |
a1f80595 | 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> | |
2180cab3 | 17 | #include <AliStack.h> |
a1f80595 | 18 | #include <iostream> |
19 | #include <TStyle.h> | |
9aa0bdc4 | 20 | #include <TArrayF.h> |
8f6ee336 | 21 | #include <TParticle.h> |
a9579262 | 22 | #include <TCanvas.h> |
23 | #include <TGraphErrors.h> | |
2b893216 | 24 | #include <TDatabasePDG.h> |
25 | #include <TParticlePDG.h> | |
26 | #include <TLegend.h> | |
27 | #include <TArrow.h> | |
28 | #include <TLatex.h> | |
29 | #include <TF1.h> | |
a1f80595 | 30 | |
9b48326f | 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 | */ | |
8ec606c2 | 41 | class DrawHits : public AliFMDInput |
a1f80595 | 42 | { |
43 | private: | |
44 | TH2D* fElossVsPMQ; // Histogram | |
2b893216 | 45 | TH1D* fEloss; |
2180cab3 | 46 | TH1D* fBetaGamma; |
2b893216 | 47 | TParticlePDG* fPdg; |
48 | const Double_t fRho; | |
49 | const Double_t fBetaGammaMip; | |
a1f80595 | 50 | public: |
bf000c32 | 51 | //__________________________________________________________________ |
2b893216 | 52 | DrawHits(const char* pdgName="pi+", |
2180cab3 | 53 | Int_t m=500, Double_t emin=1, Double_t emax=1000, |
2b893216 | 54 | Int_t n=900, Double_t tmin=1e-2, Double_t tmax=1e3) |
5632c898 | 55 | : AliFMDInput("galice.root"), |
2180cab3 | 56 | fElossVsPMQ(0), |
57 | fEloss(0), | |
58 | fBetaGamma(0), | |
5632c898 | 59 | fPdg(0), |
2180cab3 | 60 | fRho(2.33), // 2.33), |
2b893216 | 61 | fBetaGammaMip(3.4601) |
088f8e79 | 62 | { |
8f6ee336 | 63 | AddLoad(kKinematics); |
8ec606c2 | 64 | AddLoad(kHits); |
2b893216 | 65 | TDatabasePDG* pdgDB = TDatabasePDG::Instance(); |
66 | fPdg = pdgDB->GetParticle(pdgName); | |
67 | if (!fPdg) Warning("DrawHits", "Particle %s not found", pdgName); | |
68 | ||
2180cab3 | 69 | TArrayF tkine(MakeLogScale(n, tmin, tmax)); |
70 | TArrayF betag(MakeLogScale(n/10, tmin, tmax)); | |
71 | TArrayF eloss(MakeLogScale(m, emin, emax)); | |
2b893216 | 72 | TString name("elossVsPMQ"); |
73 | TString title(Form("#Delta E/#Delta x / q^{2} vs. p/m, %s", | |
74 | (pdgName ? pdgName : ""))); | |
8f6ee336 | 75 | fElossVsPMQ = new TH2D(name.Data(), title.Data(), |
76 | tkine.fN-1, tkine.fArray, | |
77 | eloss.fN-1, eloss.fArray); | |
2b893216 | 78 | fElossVsPMQ->SetXTitle("p/(mq^{2})=#beta#gamma/q^{2}"); |
79 | fElossVsPMQ->SetYTitle("#Delta E/#Delta x / q^{2} [MeV/cm]"); | |
2aeec17d | 80 | fElossVsPMQ->Sumw2(); |
2b893216 | 81 | fEloss = new TH1D("eloss", "#Delta E/#Delta x / q^{2}", |
82 | eloss.fN-1, eloss.fArray); | |
2180cab3 | 83 | fEloss->SetFillColor(kRed); |
2b893216 | 84 | fEloss->SetFillStyle(3001); |
85 | fEloss->SetXTitle("#Delta E/#Delta x / q^{2} [MeV/cm]"); | |
2aeec17d | 86 | fEloss->Sumw2(); |
2180cab3 | 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); | |
a1f80595 | 93 | } |
bf000c32 | 94 | //__________________________________________________________________ |
8f6ee336 | 95 | Bool_t ProcessHit(AliFMDHit* hit, TParticle* p) |
a1f80595 | 96 | { |
97 | if (!hit) { | |
98 | std::cout << "No hit" << std::endl; | |
99 | return kFALSE; | |
100 | } | |
088f8e79 | 101 | |
8f6ee336 | 102 | if (!p) { |
103 | std::cout << "No track" << std::endl; | |
104 | return kFALSE; | |
105 | } | |
faf80567 | 106 | // if (!p->IsPrimary()) return kTRUE; |
088f8e79 | 107 | if (hit->IsStop()) return kTRUE; |
2180cab3 | 108 | if (hit->Length() == 0) { |
109 | Warning("ProcessHit", "Hit in %s has 0 length", hit->GetName()); | |
110 | return kTRUE; | |
111 | } | |
112 | ||
2b893216 | 113 | Float_t q = hit->Q() / 3.; |
114 | Float_t m = hit->M(); | |
115 | if (m == 0 || q == 0) return kTRUE; | |
116 | ||
2180cab3 | 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 | ||
2b893216 | 140 | x /= hit->M(); |
2180cab3 | 141 | // y /= q * q; |
088f8e79 | 142 | fElossVsPMQ->Fill(x, y); |
2b893216 | 143 | fEloss->Fill(y); |
5632c898 | 144 | // fNHits++; |
a1f80595 | 145 | return kTRUE; |
146 | } | |
bf000c32 | 147 | //__________________________________________________________________ |
2180cab3 | 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 | //__________________________________________________________________ | |
a1f80595 | 178 | Bool_t Finish() |
179 | { | |
180 | gStyle->SetPalette(1); | |
2b893216 | 181 | // gStyle->SetOptTitle(0); |
182 | gStyle->SetTitleBorderSize(1); | |
183 | gStyle->SetTitleFillColor(0); | |
184 | gStyle->SetCanvasColor(0); | |
088f8e79 | 185 | gStyle->SetCanvasColor(0); |
186 | gStyle->SetCanvasBorderSize(0); | |
187 | gStyle->SetPadColor(0); | |
188 | gStyle->SetPadBorderSize(0); | |
2180cab3 | 189 | gStyle->SetOptStat(0); |
190 | TCanvas* c = new TCanvas("elossVsP", "Energy loss versus momentum", | |
191 | 1200, 800); | |
2b893216 | 192 | c->SetLogy(); |
193 | c->SetLogx(); | |
194 | ||
69893a66 | 195 | TString title(Form("%s, %d events", fElossVsPMQ->GetTitle(), fEventCount)); |
2b893216 | 196 | fElossVsPMQ->SetTitle(title.Data()); |
a1f80595 | 197 | fElossVsPMQ->SetStats(kFALSE); |
2b893216 | 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"); | |
2180cab3 | 208 | l->SetFillColor(0); |
2b893216 | 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"); | |
a9579262 | 223 | c->Modified(); |
224 | c->Update(); | |
225 | c->cd(); | |
2180cab3 | 226 | c->SaveAs("eloss_bethe.png"); |
2b893216 | 227 | |
2180cab3 | 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); | |
2b893216 | 234 | // c->SetLogx(); |
2180cab3 | 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(); | |
faf80567 | 246 | if (fEloss->GetEntries() != 0) { |
247 | c->SetLogy(); | |
248 | fEloss->Scale(1. / fEloss->GetEntries()); | |
249 | fEloss->GetXaxis()->SetRangeUser(1, 10); | |
2180cab3 | 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"); | |
faf80567 | 259 | } |
260 | ||
2180cab3 | 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); | |
2b893216 | 284 | l->AddEntry(fEloss, fEloss->GetTitle(), "lf"); |
2180cab3 | 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"); | |
2b893216 | 295 | l->Draw("same"); |
2180cab3 | 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); | |
2b893216 | 300 | |
2180cab3 | 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 | ||
a1f80595 | 334 | return kTRUE; |
335 | } | |
2b893216 | 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); | |
2180cab3 | 362 | graph->SetLineColor(kRed+1); |
363 | graph->SetLineStyle(2); | |
364 | graph->SetLineWidth(2); | |
2b893216 | 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() | |
a9579262 | 397 | { |
2b893216 | 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); | |
2180cab3 | 406 | graph->SetLineColor(kGreen+1); |
407 | graph->SetLineStyle(2); | |
408 | graph->SetLineWidth(2); | |
2b893216 | 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); | |
2180cab3 | 448 | graph->SetLineColor(kBlue+1); |
449 | graph->SetLineWidth(2); | |
450 | graph->SetLineStyle(2); | |
2b893216 | 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"); | |
69893a66 | 483 | graph->SetTitle("Mean #Delta E/#Delta x - " |
484 | "electronic only [RPP fig. 27.6]"); | |
2b893216 | 485 | graph->GetHistogram()->SetYTitle("(MeVcm^{2}/g)"); |
486 | graph->GetHistogram()->SetXTitle("#mu E_{kin} (GeV)"); | |
2180cab3 | 487 | graph->SetFillColor(0); |
488 | graph->SetLineStyle(2); | |
489 | graph->SetLineWidth(2); | |
490 | graph->SetLineColor(kMagenta+1); | |
2b893216 | 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]"); | |
2180cab3 | 523 | gre->SetFillColor(kGray); |
524 | gre->SetFillStyle(3001); // 0 /* 3002 */); | |
525 | gre->SetLineColor(kGray+1); | |
2b893216 | 526 | gre->SetLineStyle(1); |
2180cab3 | 527 | gre->SetLineWidth(2); |
2b893216 | 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 | } | |
2180cab3 | 624 | gre->Draw("c3 same"); |
2b893216 | 625 | return gre; |
626 | } | |
627 | ||
628 | /** Get the response functin @f$ f(\Delta_p/x)@f$ from Review of | |
2180cab3 | 629 | Particle Physics (fig. 27.7). It is scaled to the value at |
2b893216 | 630 | MPV. */ |
631 | TGraph* GetResp() | |
632 | { | |
633 | static TGraph* graph = 0; | |
634 | if (!graph) { | |
2180cab3 | 635 | graph = new TGraph; |
2b893216 | 636 | graph->SetName("si_resp"); |
2180cab3 | 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); | |
2b893216 | 643 | graph->SetMarkerStyle(21); |
644 | graph->SetMarkerSize(0.6); | |
2180cab3 | 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 | |
2b893216 | 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); | |
2180cab3 | 713 | #endif |
2b893216 | 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; | |
a9579262 | 749 | } |
750 | ||
8f6ee336 | 751 | ClassDef(DrawHits,0); |
a1f80595 | 752 | }; |
753 | ||
754 | //____________________________________________________________________ | |
755 | // | |
756 | // EOF | |
757 | // |