| 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 | // |
git://git.uio.no / u / mrichter / AliRoot.git / blame