| 2b893216 |
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> |
| 2180cab3 |
31 | #include <TLatex.h> |
| 2b893216 |
32 | #include <TLine.h> |
| 2180cab3 |
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 | } |
| 2b893216 |
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; |
| 2180cab3 |
76 | TList fCleanup; |
| 2b893216 |
77 | public: |
| 2b893216 |
78 | //__________________________________________________________________ |
| 2180cab3 |
79 | DrawESD(Int_t n=2000, Double_t mmin=-0.5, Double_t mmax=15.5) |
| 2b893216 |
80 | : fCorr(1) // 0.68377 / 1.1) |
| 81 | { |
| 82 | AddLoad(kESD); |
| 2180cab3 |
83 | fMult = new TH1D("mult", "#DeltaE/#DeltaE_{MIP)", n, mmin, mmax); |
| 69893a66 |
84 | fMult->Sumw2(); |
| 2180cab3 |
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(); |
| 2b893216 |
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 | //__________________________________________________________________ |
| 2180cab3 |
103 | Bool_t ProcessESD(UShort_t /* det */, |
| 104 | Char_t /* rng */, |
| 105 | UShort_t /* sec */, |
| 106 | UShort_t /* str */, |
| 107 | Float_t eta, |
| 2b893216 |
108 | Float_t mult) |
| 109 | { |
| 110 | // Cache the energy loss |
| 2180cab3 |
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 "); |
| 2b893216 |
130 | return kTRUE; |
| 131 | } |
| f48d9b11 |
132 | //__________________________________________________________________ |
| 133 | TF1* FitPeak(Int_t n, TH1D* hist, Double_t min, Double_t& max) |
| 134 | { |
| faf80567 |
135 | if (TMath::Abs(max-min) < .25) return 0; |
| 2180cab3 |
136 | std::cout << "Fit peak in range " << min << " to " << max << std::endl; |
| f48d9b11 |
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 | |
| 2180cab3 |
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 | } |
| 2b893216 |
209 | //__________________________________________________________________ |
| 210 | Bool_t Finish() |
| 211 | { |
| faf80567 |
212 | Info("Finish", "Will draw results"); |
| 2b893216 |
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 | |
| faf80567 |
221 | if (fMult->GetEntries() <= 0) return kFALSE; |
| 222 | |
| 2180cab3 |
223 | TCanvas* c = new TCanvas("c", "C", 1200, 800); |
| 224 | fCleanup.Add(c); |
| 2b893216 |
225 | c->cd(); |
| 2180cab3 |
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); |
| faf80567 |
243 | |
| 2180cab3 |
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 | |
| faf80567 |
257 | return kTRUE; |
| 258 | |
| 2180cab3 |
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 | |
| f48d9b11 |
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); |
| faf80567 |
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 | |
| 2b893216 |
494 | fMult->GetXaxis()->SetRangeUser(0, 4); |
| 2180cab3 |
495 | fMult->Fit(f, "NQR", "", x1, x3); |
| 496 | fMult->Fit(f, "NMER", "E1", x1, x3); |
| 2b893216 |
497 | |
| 2b893216 |
498 | l1->SetLineColor(3); |
| 499 | l1->SetLineWidth(2); |
| 500 | l1->SetRange(0, 4); |
| 501 | l1->Draw("same"); |
| faf80567 |
502 | if (l2) { |
| 503 | l2->SetLineColor(4); |
| 504 | l2->SetLineWidth(2); |
| 505 | l2->SetRange(0, 4); |
| 506 | l2->Draw("same"); |
| 507 | } |
| 2b893216 |
508 | f->SetLineWidth(2); |
| 509 | f->SetRange(0, 4); |
| 510 | f->Draw("same"); |
| 511 | |
| 2180cab3 |
512 | fCleanup.Add(l1); |
| 513 | if (l2) fCleanup.Add(l2); |
| 514 | fCleanup.Add(f); |
| 2b893216 |
515 | |
| 516 | |
| 2180cab3 |
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"); |
| faf80567 |
532 | // c->SetLogy(); |
| 2b893216 |
533 | fMult->DrawClone("axis"); |
| 534 | f->Draw("same"); |
| f48d9b11 |
535 | Double_t* p1 = f->GetParameters(); |
| 536 | Double_t* p2 = &(p1[3]); |
| 2b893216 |
537 | TF1* ll1 = new TF1("ll1", "landau", 0, 4); |
| f48d9b11 |
538 | ll1->SetParameters(p1); |
| 2b893216 |
539 | ll1->SetLineColor(3); |
| 540 | ll1->Draw("same"); |
| 541 | TF1* ll2 = new TF1("ll2", "landau", 0, 4); |
| f48d9b11 |
542 | ll2->SetParameters(p2); |
| 2b893216 |
543 | ll2->SetLineColor(4); |
| 544 | ll2->Draw("same"); |
| 545 | |
| 2180cab3 |
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; |
| 2b893216 |
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"); |
| f48d9b11 |
556 | TLine* c3 = new TLine(xc3, y1, xc3, y2); |
| 557 | c3->Draw("same"); |
| 2b893216 |
558 | |
| 2180cab3 |
559 | TLegend* l = new TLegend(0.6, 0.6, .89, .89); |
| 2b893216 |
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"); |
| 2b893216 |
565 | |
| 2180cab3 |
566 | c->Modified(); |
| 567 | c->Update(); |
| 568 | c->cd(); |
| 569 | } |
| f48d9b11 |
570 | |
| 2180cab3 |
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; |
| 2b893216 |
747 | } |
| 748 | |
| 749 | ClassDef(DrawESD,0); |
| 750 | |
| 751 | }; |
| 752 | |
| 753 | //____________________________________________________________________ |
| 754 | // |
| 755 | // EOF |
| 756 | // |
git://git.uio.no / u / mrichter / AliRoot.git / blame