]>
Commit | Line | Data |
---|---|---|
3cd4fae8 | 1 | #include "AliLandauGaus.h" |
2 | #include <TSystem.h> | |
3 | #include <TString.h> | |
4 | #include <TList.h> | |
5 | #include <TFile.h> | |
6 | #include <TMath.h> | |
7 | #include <TF1.h> | |
8 | #include <TH1.h> | |
9 | #include <TH2.h> | |
10 | #include <TFitResult.h> | |
11 | #include <TNtuple.h> | |
12 | #include <TAxis.h> | |
13 | #include <TMath.h> | |
14 | #include <TCanvas.h> | |
15 | #include <TStyle.h> | |
16 | #include <TLine.h> | |
17 | #include <TLatex.h> | |
18 | #include <TPaveStats.h> | |
19 | ||
20 | //==================================================================== | |
21 | /** | |
22 | * Base class for tests | |
23 | * | |
24 | */ | |
25 | struct TestF | |
26 | { | |
27 | //__________________________________________________________________ | |
28 | /** | |
29 | * Base class for individual tests | |
30 | */ | |
31 | struct Test { | |
32 | /** Source distribution */ | |
33 | TF1* src; | |
34 | /** | |
35 | * Constructor | |
36 | */ | |
37 | Test() : src(0) {} | |
38 | /** | |
39 | * Copy constructor | |
40 | */ | |
41 | Test(const Test& o) | |
42 | : src(static_cast<TF1*>(o.src ? o.src->Clone() : 0)) {} | |
43 | /** | |
44 | * Destructor | |
45 | */ | |
46 | virtual ~Test() { delete src; } | |
47 | /** | |
48 | * Assignment operator | |
49 | */ | |
50 | Test& operator=(const Test& o) { | |
51 | if (&o == this) return *this; | |
52 | if (o.src) src = static_cast<TF1*>(o.src->Clone()); | |
53 | return *this; | |
54 | } | |
55 | /** | |
56 | * Run the test | |
57 | * | |
58 | * @param c Constant | |
59 | * @param delta Most probable value | |
60 | * @param xi Landau width | |
61 | * @param sigma Gaussian smear | |
62 | * @param n Maximum number of particles | |
63 | * @param a Particle weights | |
64 | * @param xMin Least energy loss | |
65 | * @param xMax Largest energy loss | |
66 | */ | |
67 | void Run(const Double_t c, | |
68 | const Double_t delta, | |
69 | const Double_t xi, | |
70 | const Double_t sigma, | |
71 | const Int_t n, | |
72 | const Double_t* a, | |
73 | const Double_t xMin, | |
74 | const Double_t xMax) { | |
75 | printf("%8.4f | %8.4f | %8.4f | ", delta, xi, sigma); | |
76 | fflush(stdout); | |
77 | src = AliLandauGaus::MakeFn(c,delta,xi,sigma,0,n,&(a[1]),xMin,xMax); | |
78 | src->SetName("src"); | |
79 | src->SetLineColor(kBlack); | |
80 | src->SetNpx(10000); | |
81 | // If the constant is 1, then this is a fake distribtion and we | |
82 | // need to scale it. | |
83 | if (c >= 1) { | |
84 | Double_t max = src->GetMaximum(); | |
85 | src->SetParameter(0, c/max); | |
86 | } | |
87 | ||
88 | // src->Print(); | |
89 | DoRun(); | |
90 | } | |
91 | /** | |
92 | * @return Source parameter @f$ C@f$ | |
93 | */ | |
94 | Double_t C() const { return src->GetParameter(AliLandauGaus::kC); } | |
95 | /** | |
96 | * @return Source parameter @f$\Delta_p@f$ | |
97 | */ | |
98 | Double_t Delta() const { return src->GetParameter(AliLandauGaus::kDelta); } | |
99 | /** | |
100 | * @return Source parameter @f$\xi@f$ | |
101 | */ | |
102 | Double_t Xi() const { return src->GetParameter(AliLandauGaus::kXi); } | |
103 | /** | |
104 | * @return Source parameter @f$ \sigma@f$ | |
105 | */ | |
106 | Double_t Sigma() const { return src->GetParameter(AliLandauGaus::kSigma); } | |
107 | /** | |
108 | * @return Source parameter @f$ \sigma_n@f$ | |
109 | */ | |
110 | Double_t SigmaN() const { return src->GetParameter(AliLandauGaus::kSigmaN);} | |
111 | /** | |
112 | * @return Source parameter @f$ N@f$ | |
113 | */ | |
114 | Double_t N() const { return src->GetParameter(AliLandauGaus::kN); } | |
115 | /** | |
116 | * @return Source parameter @f$ a_i@f$ | |
117 | */ | |
118 | Double_t A(Int_t i) const { | |
119 | return (i <= 1 ? 1 : src->GetParameter(AliLandauGaus::kN+(i-1))); } | |
120 | /** | |
121 | * @return Least energy loss | |
122 | */ | |
123 | Double_t XMin() const { return src->GetXmin(); } | |
124 | /** | |
125 | * @return Largest energy loss | |
126 | */ | |
127 | Double_t XMax() const { return src->GetXmax(); } | |
128 | /** | |
129 | * Generate single component function for @a i particles | |
130 | * | |
131 | * @param c Constant | |
132 | * @param delta @f$\Delta_p@f$ | |
133 | * @param xi @f$\xi@f$ | |
134 | * @param sigma @f$\sigma@f$ | |
135 | * @param sigmaN @f$\sigma_n@f$ | |
136 | * @param i Number of particles | |
137 | * | |
138 | * @return Pointer to function object | |
139 | */ | |
140 | TF1* Fi(Double_t c, Double_t delta, Double_t xi, Double_t sigma, | |
141 | Double_t sigmaN, Int_t i) const { | |
142 | return AliLandauGaus::MakeFi(c,delta,xi,sigma,sigmaN,i,XMin(),XMax()); | |
143 | } | |
144 | /** | |
145 | * Generate single component function for @a i particles using | |
146 | * source parameters | |
147 | * | |
148 | * @param i Number of particles | |
149 | * | |
150 | * @return Pointer to function object | |
151 | */ | |
152 | TF1* Fi(Int_t i) { | |
153 | return Fi(C()*A(i),Delta(),Xi(),Sigma(),SigmaN(),i); | |
154 | } | |
155 | ||
156 | /** | |
157 | * @{ | |
158 | * @name Overload-able interface | |
159 | */ | |
160 | /** | |
161 | * Run the test | |
162 | */ | |
163 | virtual void DoRun() = 0; | |
164 | /** | |
165 | * Write results to disk | |
166 | * | |
167 | * @param d Directory to write to | |
168 | */ | |
169 | virtual void WriteOut(TDirectory* d) = 0; | |
170 | /** | |
171 | * Draw results in a pad | |
172 | * | |
173 | * @param p Pad to draw in | |
174 | */ | |
175 | virtual void DrawInPad(TVirtualPad* p) = 0; | |
176 | /* @} */ | |
177 | }; | |
178 | //__________________________________________________________________ | |
179 | /** | |
180 | * @{ | |
181 | * @name Utilities | |
182 | */ | |
183 | /** | |
184 | * Draw a vertical line | |
185 | * | |
186 | * @param i Line style | |
187 | * @param col Line color | |
188 | * @param x X value | |
189 | * @param low Least Y value | |
190 | * @param high Largest Y value | |
191 | * | |
192 | * @return The line | |
193 | */ | |
194 | static TLine* | |
195 | VerticalLine(Int_t i, Color_t col, Double_t x, Double_t low, Double_t high) | |
196 | { | |
197 | TLine* l = new TLine(x,low,x,high); | |
198 | l->SetLineStyle(i); | |
199 | l->SetLineColor(col); | |
200 | return l; | |
201 | } | |
202 | //__________________________________________________________________ | |
203 | virtual ~TestF() {} | |
204 | /** | |
205 | * Make a test object - overload | |
206 | * | |
207 | * @return The test object | |
208 | */ | |
209 | virtual Test* MakeTest() = 0; | |
210 | /** | |
211 | * Do one scan | |
212 | * | |
213 | * @param c Constant | |
214 | * @param delta Most probable value @f$\Delta_p@f$ | |
215 | * @param xi Landau width @f$\xi@f$ | |
216 | * @param sigma Gaussian spread @f$\sigma@f$ | |
217 | * @param n Maximum number of particles | |
218 | * @param a Particle weights | |
219 | * @param xMin Least @f$\Delta @f$ | |
220 | * @param xMax Largest @f$\Delta@f$ | |
221 | * | |
222 | * @return The test object | |
223 | */ | |
224 | Test* DoOne(const Double_t c, | |
225 | const Double_t delta, | |
226 | const Double_t xi, | |
227 | const Double_t sigma, | |
228 | const Int_t n, | |
229 | const Double_t* a, | |
230 | const Double_t xMin, | |
231 | const Double_t xMax, | |
232 | TVirtualPad* pad=0) | |
233 | { | |
234 | Test* t = MakeTest(); | |
235 | t->Run(c, delta, xi, sigma, n, a, xMin, xMax); | |
236 | ||
237 | TString name(Form("sigma%04d_xi%04d", Int_t(100*sigma), Int_t(100*xi))); | |
238 | TString title(Form("Sigma=%f xi=%f", sigma, xi)); | |
239 | if (!pad) { | |
240 | TCanvas* can = new TCanvas(name, title, 1200, 900); | |
241 | pad = can; | |
242 | } | |
243 | else | |
244 | pad->Clear(); | |
245 | t->DrawInPad(pad); | |
246 | ||
247 | return t; | |
248 | } | |
249 | /** | |
250 | * Do a unit scan | |
251 | * | |
252 | * @param sigma @f$\sigma@f$ of the Gaussian | |
253 | * @param xi @f$\xi@f$ width of Landau | |
254 | * @param n Number of particles | |
255 | * @param pad Possibly existing canvas to plot in | |
256 | * | |
257 | * @return Scan object | |
258 | */ | |
259 | Test* DoUnit(const Double_t sigma, | |
260 | const Double_t xi=1, | |
261 | const Int_t n=10, | |
262 | TVirtualPad* pad=0) | |
263 | { | |
264 | const Double_t c = 1; | |
265 | const Double_t delta = 10; | |
266 | Double_t idelta = delta; | |
267 | Double_t ixi = xi; | |
268 | Double_t isigma = sigma; | |
269 | AliLandauGaus::IPars(n+1,idelta,ixi, isigma); | |
270 | Double_t up = idelta + ixi + isigma; // 8*n | |
271 | Double_t a[n]; | |
272 | a[0] = 1; | |
273 | for (Int_t i = 1; i < n; i++) | |
274 | a[i] = TMath::Exp(0.5 - TMath::Sqrt(2)*(i+1)); | |
275 | ||
276 | return DoOne(c,delta,xi,sigma,n,a,-3,up, pad); | |
277 | } | |
278 | /** | |
279 | * Do a realistic scan | |
280 | * | |
281 | * @param pad Possibly existing canvas to plot in | |
282 | * | |
283 | * @return Scan object | |
284 | */ | |
285 | Test* DoRealistic(TVirtualPad* pad=0) | |
286 | { | |
287 | Double_t c = 0.4070; | |
288 | Double_t delta = 0.5651; | |
289 | Double_t xi = 0.0643; | |
290 | Double_t sigma = 0.0611; | |
291 | Int_t n = 5; | |
292 | Double_t a[] = {1, 0.1179, 0.0333, 0.0068, 0.0012 }; | |
293 | ||
294 | return DoOne(c, delta, xi, sigma, n, a, 0.02, n, pad); | |
295 | } | |
296 | /** | |
297 | * Prefix on output | |
298 | * | |
299 | * | |
300 | * @return | |
301 | */ | |
302 | virtual const char* Prefix() const = 0; | |
303 | /** | |
304 | * Make an output ntuple | |
305 | * | |
306 | * @return Possibly newly allocated NTuple | |
307 | */ | |
308 | virtual TNtuple* MakeNTuple() { return 0; } | |
309 | /** | |
310 | * Setup for scanning | |
311 | * | |
312 | * @param mode Scanning mode | |
313 | * @param file Output file | |
314 | * @param canvas Canvas | |
315 | * @param nt Possible Ntuple | |
316 | */ | |
317 | virtual void SetupScan(UShort_t mode, | |
318 | TFile*& file, | |
319 | TCanvas*& canvas, | |
320 | TNtuple*& nt) | |
321 | { | |
322 | TString base = Form("%s%s", Prefix(), | |
323 | (mode == 0 ? "SigmaXi" : | |
324 | mode == 1 ? "Sigma" : "Xi")); | |
325 | file = TFile::Open(Form("%s.root", base.Data()), "RECREATE"); | |
326 | nt = MakeNTuple(); | |
327 | canvas = new TCanvas(base, Form("%s.pdf", base.Data()), | |
328 | TMath::Sqrt(2)*900, 900); | |
329 | TString title("pdf Landscape Title:Scan"); | |
330 | gSystem->RedirectOutput("/dev/null"); | |
331 | canvas->Print(Form("%s[", canvas->GetTitle()), title); | |
332 | gSystem->RedirectOutput(0); | |
333 | } | |
334 | /** | |
335 | * Print canvas to output PDF | |
336 | * | |
337 | * @param c Canvas | |
338 | * @param title Title of page | |
339 | * | |
340 | * @return The full title of the page | |
341 | */ | |
342 | void PrintCanvas(TCanvas* c, Double_t xi, Double_t sigma) { | |
343 | TString tit = "pdf Landscape Title:"; | |
344 | Bool_t cls = false; | |
345 | if (xi > 0 && sigma > 0) { | |
346 | tit.Append(Form("xi=%8.4f sigma=%8.4f", xi, sigma)); | |
347 | } | |
348 | else { | |
349 | tit.Append("Results"); | |
350 | cls = true; | |
351 | } | |
352 | gSystem->RedirectOutput("/dev/null"); | |
353 | c->Print(c->GetTitle(), tit); | |
354 | if (cls) { | |
355 | c->Clear(); | |
356 | c->Print(Form("%s]", c->GetTitle()), tit); | |
357 | } | |
358 | gSystem->RedirectOutput(0); | |
359 | if (cls) | |
360 | ::Info("", "Plots saved in %s", c->GetTitle()); | |
361 | } | |
362 | /** | |
363 | * Executed before loop | |
364 | * | |
365 | * @param mode Mode | |
366 | * @param nVal Number of values | |
367 | */ | |
368 | virtual void PreLoop(UShort_t mode, Int_t nVal) = 0; | |
369 | /** | |
370 | * Process a single step | |
371 | * | |
372 | * @param mode Mode | |
373 | * @param i First index | |
374 | * @param j Second index (if any) | |
375 | * @param n Number of values | |
376 | * @param test Test result | |
377 | */ | |
378 | virtual void Step(UShort_t mode, | |
379 | Int_t i, | |
380 | Int_t j, | |
381 | Int_t n, | |
382 | Test* test, | |
383 | TNtuple* nt) = 0; | |
384 | /** | |
385 | * Called after looping | |
386 | * | |
387 | * @param mode Execution mode | |
388 | * @param can Canvas | |
389 | * @param nt NTuple | |
390 | * @param out Output file | |
391 | */ | |
392 | virtual void PostLoop(UShort_t mode, | |
393 | TCanvas* can, | |
394 | TNtuple* nt, | |
395 | TFile* out) = 0; | |
396 | /** | |
397 | * Scan over a parameter (@f$\xi@f$ or @f$\sigma@f$) | |
398 | * | |
399 | * @param sigma If true, scan @f$\sigma@f$, otherwise @f$\xi@f$ | |
400 | * @param n Number of values | |
401 | * @param values Values of the parameters | |
402 | * @param maxN Number of particles | |
403 | */ | |
404 | void ScanOne(Bool_t scanSigma, | |
405 | Int_t n, | |
406 | const Double_t* values, | |
407 | Int_t maxN=10) | |
408 | { | |
409 | UShort_t mod = scanSigma ? 1 : 2; | |
410 | Scan(mod, n, values, maxN); | |
411 | } | |
412 | /** | |
413 | * Scan over a parameter (@f$\xi@f$ or @f$\sigma@f$) | |
414 | * | |
415 | * @param sigma If true, scan @f$\sigma@f$, otherwise @f$\xi@f$ | |
416 | * @param values Values of the parameters | |
417 | * @param maxN Number of particles | |
418 | */ | |
419 | void ScanOne(Bool_t scanSigma, | |
420 | const TArrayD& values, | |
421 | Int_t maxN) | |
422 | { | |
423 | ScanOne(scanSigma, values.GetSize(), values.GetArray(), maxN); | |
424 | } | |
425 | /** | |
426 | * Scan over both parameter (@f$\xi@f$ and @f$\sigma@f$) | |
427 | * | |
428 | * @param n Number of parameters | |
429 | * @param values Values of the parameters | |
430 | * @param maxN Number of particles | |
431 | */ | |
432 | void ScanTwo(Int_t n, | |
433 | const Double_t* values, | |
434 | Int_t maxN=10) | |
435 | { | |
436 | Scan(0, n, values, maxN); | |
437 | } | |
438 | /** | |
439 | * Scan over both parameter (@f$\xi@f$ and @f$\sigma@f$) | |
440 | * | |
441 | * @param values Values of the parameters | |
442 | * @param maxN Number of particles | |
443 | */ | |
444 | void ScanTwo(const TArrayD& values, Int_t maxN) | |
445 | { | |
446 | ScanTwo(values.GetSize(), values.GetArray(), maxN); | |
447 | } | |
448 | /** | |
449 | * Do the scan | |
450 | * | |
451 | * @param mode Mode of operation (0: both, 1: @f$\sigma@f$, or @f$\xi@f$) | |
452 | * @param n Number of values | |
453 | * @param values Values | |
454 | * @param maxN Maximum number of particles | |
455 | */ | |
456 | void Scan(UShort_t mode, | |
457 | Int_t n, | |
458 | const Double_t* values, | |
459 | Int_t maxN) | |
460 | { | |
461 | TFile* out = 0; | |
462 | TNtuple* nt = 0; | |
463 | TCanvas* can = 0; | |
464 | SetupScan(mode, out, can, nt); | |
465 | PreLoop(mode, n); | |
466 | ||
467 | Test* t = DoRealistic(can); | |
468 | t->WriteOut(out->mkdir("realistic")); | |
469 | can->Write(); | |
470 | PrintCanvas(can, t->Xi(), t->Sigma()); | |
471 | ||
472 | if (mode == 0) LoopTwo(n, values, maxN, can, nt, out); | |
473 | else LoopOne(mode, n, values, maxN, can, nt, out); | |
474 | ||
475 | PostLoop(mode, can, nt, out); | |
476 | out->Write(); | |
477 | out->Close(); | |
478 | ||
479 | } | |
480 | /** | |
481 | * Loop over one of @f$\xi@f$ or @f$\sigma@f$ | |
482 | * | |
483 | * @param mode Mode (1: @f$\sigma@f$, otherwise @f$\xi@f$) | |
484 | * @param n Number of values | |
485 | * @param values Values | |
486 | * @param maxN Maximum number of particles | |
487 | * @param can Canvas to draw in | |
488 | * @param nt NTuple to fill | |
489 | * @param out Output directory | |
490 | */ | |
491 | void LoopOne(UShort_t mode, | |
492 | Int_t n, | |
493 | const Double_t* values, | |
494 | Int_t maxN, | |
495 | TCanvas* can, | |
496 | TNtuple* nt, | |
497 | TFile* out) { | |
498 | Bool_t scanSigma = mode == 1; | |
499 | const char* var = (scanSigma ? "sigma" : "xi"); | |
500 | ||
501 | for (Int_t i = 0; i < n; i++) { | |
502 | Double_t v = values[i]; | |
503 | Double_t sigma = (scanSigma ? v : 1); | |
504 | Double_t xi = (!scanSigma ? v : 1); | |
505 | ||
506 | StepIt(mode, xi, sigma, maxN, 0, i, n, can, nt, out, var, v); | |
507 | } | |
508 | } | |
509 | /** | |
510 | * Loop over both @f$\xi@f$ and @f$\sigma@f$ | |
511 | * | |
512 | * @param n Number of values | |
513 | * @param values Values | |
514 | * @param maxN Maximum number of particles | |
515 | * @param can Canvas to draw in | |
516 | * @param nt NTuple to fill | |
517 | * @param out Output directory | |
518 | */ | |
519 | void LoopTwo(Int_t n, | |
520 | const Double_t* values, | |
521 | Int_t maxN, | |
522 | TCanvas* can, | |
523 | TNtuple* nt, | |
524 | TFile* out) | |
525 | { | |
526 | UShort_t mode = 0; | |
527 | for (Int_t i = 0; i < n; i++) { | |
528 | Double_t xi = values[i]; | |
529 | TDirectory* d = out->mkdir(Form("xi%04d", Int_t(100*xi))); | |
530 | for (Int_t j = 0; j < n; j++) { | |
531 | Double_t sigma = values[j]; | |
532 | StepIt(mode, xi, sigma, maxN, i, j, n, can, nt, d, "sigma", sigma); | |
533 | } | |
534 | } | |
535 | } | |
536 | /** | |
537 | * | |
538 | * | |
539 | * @param mode Mode | |
540 | * @param xi @f$\xi@f$ | |
541 | * @param sigma @f$\sigma@f$ | |
542 | * @param maxN At most this many particles | |
543 | * @param i First index | |
544 | * @param j Last index | |
545 | * @param n Number of values | |
546 | * @param can Canvas | |
547 | * @param nt NTuple | |
548 | * @param p Parent directory | |
549 | * @param pre Prefix on directory | |
550 | * @param v Current value | |
551 | */ | |
552 | void StepIt(UShort_t mode, | |
553 | Double_t xi, | |
554 | Double_t sigma, | |
555 | Int_t maxN, | |
556 | Int_t i, | |
557 | Int_t j, | |
558 | Int_t n, | |
559 | TCanvas* can, | |
560 | TNtuple* nt, | |
561 | TDirectory* p, | |
562 | const char* pre, | |
563 | Double_t v) | |
564 | { | |
565 | TDirectory* d = p->mkdir(Form("%s%04d", pre, Int_t(100*v))); | |
566 | Test* t = DoUnit(sigma, xi, maxN, can); | |
567 | t->WriteOut(d); | |
568 | can->Write(); | |
569 | PrintCanvas(can, xi, sigma); | |
570 | ||
571 | Step(mode, i, j, n, t, nt); | |
572 | } | |
573 | }; | |
574 | ||
575 | //==================================================================== | |
576 | #ifdef TEST_SHIFT | |
577 | #include "WrappedGraph.C" | |
578 | #include <TGraphErrors.h> | |
579 | #include <TGraph2DErrors.h> | |
580 | #include <TMultiGraph.h> | |
581 | #include <Math/RootFinder.h> | |
582 | ||
583 | struct TestShift : public TestF | |
584 | { | |
585 | /** | |
586 | * Our fit function | |
587 | * | |
588 | * @f[ | |
589 | * \Delta_{i,p} - \Delta_{i,r} = | |
590 | * f(i;c,p,\xi,\sigma) = | |
591 | * \frac{c u \sigma}{(1+1/i)^{p u^{3/2}} | |
592 | * @f] | |
593 | * | |
594 | * where @f$ u=\sigma/\xi@f$ | |
595 | * | |
596 | * @param xp Independent variables | |
597 | * @param pp Parameters | |
598 | * | |
599 | * @return @f$ f(i;c,p,\xi,\sigma)@f$ | |
600 | */ | |
601 | static Double_t fitFunc(Double_t* xp, Double_t* pp) { | |
602 | Double_t x = *xp; | |
603 | Double_t c = pp[0]; | |
604 | Double_t p = pp[1]; | |
605 | Double_t xi = pp[2]; | |
606 | Double_t sigma = pp[3]; | |
607 | Double_t u = sigma / xi; | |
608 | Double_t a = c * u * sigma; | |
609 | Double_t b = p * u * TMath::Sqrt(u); | |
610 | ||
611 | return a / TMath::Power(1+1./x, b); | |
612 | } | |
613 | /** | |
614 | * Scale a graph | |
615 | * | |
616 | * @param g Graph to scale | |
617 | * @param scale Scale | |
618 | */ | |
619 | static void ScaleGraph(TGraph* g, Double_t scale) | |
620 | { | |
621 | // Info("", "Scaling graph %p with %f", g, scale); | |
622 | if (scale == 1) return; | |
623 | for (Int_t i = 0; i < g->GetN(); i++) | |
624 | g->SetPoint(i,g->GetX()[i], g->GetY()[i] * scale); | |
625 | } | |
626 | //__________________________________________________________________ | |
627 | struct Test : public TestF::Test | |
628 | { | |
629 | TGraph* zeros; // Zeros of the derivatives | |
630 | TGraph* diff; // Difference from calculated Delta to roots | |
631 | TList comps; // Components | |
632 | TList derivs; // Derivatives | |
633 | TF1* fit; // Fit to difference | |
634 | Test() : TestF::Test(), zeros(0), diff(0), fit(0) | |
635 | { | |
636 | comps.SetName("components"); | |
637 | derivs.SetName("derivatives"); | |
638 | ||
639 | } | |
640 | /** | |
641 | * Run the actual test | |
642 | */ | |
643 | void DoRun() | |
644 | { | |
645 | AliLandauGaus::EnableSigmaShift(false); // Disable sigma-shift | |
646 | const Int_t n = N(); | |
647 | zeros = new TGraph(n); | |
648 | zeros->SetName("zeros"); | |
649 | zeros->SetTitle("Zeros of derivatives"); | |
650 | zeros->SetMarkerStyle(25); | |
651 | zeros->SetMarkerSize(1.2); | |
652 | zeros->SetMarkerColor(kRed+1); | |
653 | ||
654 | diff = new TGraph(n); | |
655 | diff->SetName("diff"); | |
656 | diff->SetTitle("#delta#Delta_{p}"); | |
657 | diff->SetMarkerStyle(20); | |
658 | diff->SetMarkerColor(kBlue+1); | |
659 | ||
660 | // Make the components | |
661 | Int_t j = 0; | |
662 | for (Int_t i = 1; i <= n; i++, j++) | |
663 | if (!MakeComp(i)) break; | |
664 | zeros->Set(j); | |
665 | diff->Set(j); | |
666 | ||
667 | // Now fit | |
668 | fit = new TF1("fit", &fitFunc, 1, j, 4); | |
669 | fit->SetParNames("c", "p", "#xi", "#sigma"); | |
670 | fit->SetParameter(0, .5); | |
671 | fit->SetParameter(1, .5); | |
672 | fit->FixParameter(2, Xi()); | |
673 | fit->FixParameter(3, Sigma()); | |
674 | ||
675 | // fit->SetParNames("c", "#xi", "#sigma"); | |
676 | // fit->FixParameter(1, xi); | |
677 | ||
678 | ||
679 | Int_t ifit = 4; | |
680 | do { | |
681 | gSystem->RedirectOutput("/dev/null"); | |
682 | diff->Fit(fit, Form("MR%s", ifit == 0 ? "Q" : "Q")); | |
683 | gSystem->RedirectOutput(0); | |
684 | ifit--; | |
685 | } while (ifit >= 0); | |
686 | Printf("%8.4f | %8.4f %9.5f | %8.4f %9.5f", | |
687 | fit->GetChisquare() / fit->GetNDF(), | |
688 | fit->GetParameter(0), fit->GetParError(0), | |
689 | fit->GetParameter(1), fit->GetParError(1)); | |
690 | ||
691 | } | |
692 | /** | |
693 | * Make a single component | |
694 | * | |
695 | * @param i Number of particles | |
696 | */ | |
697 | Bool_t MakeComp(const Int_t i) { | |
698 | TF1* fi = Fi(i); | |
699 | fi->SetNpx(src->GetNpx()); | |
700 | // fi->Print(); | |
701 | ||
702 | TGraph* dF = new TGraph(fi, "d"); | |
703 | dF->SetName(Form("d%s", fi->GetName())); | |
704 | dF->SetTitle(Form("d%s/d#Delta", fi->GetTitle())); | |
705 | dF->SetMarkerStyle(1); // A dot | |
706 | // dF->SetMarkerStyle(20+i-1); | |
707 | dF->SetMarkerColor(fi->GetLineColor()); | |
708 | dF->SetLineColor(fi->GetLineColor()); | |
709 | ||
710 | Double_t max = TMath::MaxElement(dF->GetN(), dF->GetY()); | |
711 | if (max < 1e-6) return false; | |
712 | ||
713 | ScaleGraph(dF, 1./max); | |
714 | Double_t delta = Delta(); | |
715 | Double_t xi = Xi(); | |
716 | Double_t deltaI = i * (delta + xi * TMath::Log(i)); | |
717 | Double_t xiI = i * xi; | |
718 | Double_t xR = FindRoot(dF, deltaI-2*xiI, deltaI+2*xiI); | |
719 | ||
720 | // Printf("Component %2d: Exp: %8.4f Got: %8.4f Diff: %9.5f", | |
721 | // i, deltaI, xR, (xR-deltaI)); | |
722 | // printf("."); | |
723 | comps.Add(fi); | |
724 | derivs.Add(dF); | |
725 | zeros->SetPoint(i-1, xR, dF->Eval(xR)); | |
726 | diff->SetPoint(i-1, i, xR - deltaI); | |
727 | ||
728 | return true; | |
729 | } | |
730 | /** | |
731 | * Find the root of a graph | |
732 | * | |
733 | * @param g Graph | |
734 | * @param xMin Least value of scan range | |
735 | * @param xMax Largest value of scan range | |
736 | * | |
737 | * @return Zero of the graph | |
738 | */ | |
739 | Double_t FindRoot(TGraph* g, Double_t xMin, Double_t xMax) { | |
740 | WrappedGraph* wG = new WrappedGraph(g); | |
741 | ROOT::Math::RootFinder rfb(ROOT::Math::RootFinder::kBRENT); | |
742 | rfb.SetFunction(*wG, xMin, xMax); | |
743 | rfb.Solve(); | |
744 | Double_t xR = rfb.Root(); | |
745 | return xR; | |
746 | } | |
747 | /** | |
748 | * Draw the results in a pad | |
749 | * | |
750 | * @param body Pad to draw in | |
751 | */ | |
752 | void DrawInPad(TVirtualPad* body) { | |
753 | body->SetTopMargin(0.01); | |
754 | body->SetRightMargin(0.01); | |
755 | body->Divide(2,1); | |
756 | ||
757 | TVirtualPad* p = body->cd(1); | |
758 | p->SetTopMargin(0.01); | |
759 | p->SetRightMargin(0.01); | |
760 | p->Divide(1,2,0,0); | |
761 | ||
762 | TVirtualPad* q = p->cd(1); | |
763 | q->SetRightMargin(0.01); | |
764 | q->SetGridx(); | |
765 | q->SetLogy(); | |
766 | src->Draw(); | |
767 | Double_t fmax = src->GetHistogram()->GetMaximum(); | |
768 | Double_t fmin = src->GetHistogram()->GetMinimum(); | |
769 | fmin = TMath::Max(fmin, fmax/1e8); | |
770 | src->GetHistogram()->SetMinimum(fmin); | |
771 | TLatex* ltx = new TLatex(q->GetLeftMargin()+.02, | |
772 | 1-q->GetTopMargin()-.02, | |
773 | Form("#Delta_{p}=%f, #xi=%f, #sigma=%f", | |
774 | Delta(), Xi(), Sigma())); | |
775 | // Printf("%s: Max=%f scale=%f", src->GetName(), fmax, 1/fmax); | |
776 | ltx->SetTextFont(42); | |
777 | ltx->SetNDC(); | |
778 | ltx->SetTextAlign(13); | |
779 | ltx->Draw(); | |
780 | ||
781 | q = p->cd(2); | |
782 | q->SetRightMargin(0.01); | |
783 | q->SetGridx(); | |
784 | ||
785 | Double_t x = 1-q->GetRightMargin()-.02; | |
786 | Double_t y = 1-q->GetTopMargin()-.1; | |
787 | for (Int_t i = 1; i <= comps.GetEntries(); i++) { | |
788 | p->cd(1); | |
789 | TF1* fi = static_cast<TF1*>(comps.At(i-1)); | |
790 | fi->Draw("same"); | |
791 | ||
792 | Double_t deltaR = zeros->GetX()[i-1]; | |
793 | Double_t deltaI = i * (Delta() + Xi() * TMath::Log(i)); | |
794 | ||
795 | Color_t col = fi->GetLineColor(); | |
796 | VerticalLine(1,col,deltaI,fmin,fmax)->Draw(); | |
797 | VerticalLine(3,col,deltaR,fmin,fmax)->Draw(); | |
798 | ||
799 | ltx = new TLatex(x,y,Form("#Delta_{p,%2d}=%8.4f %8.4f", | |
800 | i,deltaI,deltaR)); | |
801 | ltx->SetTextAlign(33); | |
802 | ltx->SetTextFont(42); | |
803 | ltx->SetNDC(); | |
804 | ltx->Draw(); | |
805 | y -= ltx->GetTextSize() + .01; | |
806 | ||
807 | p->cd(2); | |
808 | ||
809 | TGraph* df = static_cast<TGraph*>(derivs.At(i-1)); | |
810 | df->Draw(i == 1 ? "alp" : "lp"); | |
811 | df->GetHistogram()->GetXaxis()->SetRangeUser(XMin(), XMax()); | |
812 | VerticalLine(1,col,deltaI,-.6,1.1)->Draw(); | |
813 | VerticalLine(3,col,deltaR,-.6,1.1)->Draw(); | |
814 | } | |
815 | p->cd(2); | |
816 | zeros->Draw("p"); | |
817 | ||
818 | gStyle->SetOptFit(111111); | |
819 | gStyle->SetStatY(.6); | |
820 | p = body->cd(2); | |
821 | p->SetRightMargin(0.02); | |
822 | p->SetTopMargin(0.02); | |
823 | diff->Draw("alp"); | |
824 | diff->GetHistogram()->GetListOfFunctions()->Add(fit); | |
825 | diff->GetHistogram()->SetXTitle("N_{particle}"); | |
826 | diff->GetHistogram()->SetYTitle("#Delta_{r}-#Delta_{p}"); | |
827 | p->Clear(); | |
828 | diff->Draw("alp"); | |
829 | ||
830 | body->Modified(); | |
831 | body->Update(); | |
832 | body->cd(); | |
833 | } | |
834 | /** | |
835 | * Write results to disk | |
836 | * | |
837 | * @param dir Directory | |
838 | */ | |
839 | void WriteOut(TDirectory* dir) { | |
840 | dir->cd(); | |
841 | src->Write(); | |
842 | zeros->Write(); | |
843 | diff->Write(); | |
844 | comps.Write(comps.GetName(), TObject::kSingleKey); | |
845 | derivs.Write(derivs.GetName(), TObject::kSingleKey); | |
846 | fit->Write(); | |
847 | } | |
848 | }; | |
849 | //__________________________________________________________________ | |
850 | TGraphErrors* cs; | |
851 | TGraphErrors* ps; | |
852 | TGraph2DErrors* c2s; | |
853 | TGraph2DErrors* p2s; | |
854 | TMultiGraph* mcs; | |
855 | TMultiGraph* mps; | |
856 | ||
857 | /** | |
858 | * Make a test object - overload | |
859 | * | |
860 | * @return The test object | |
861 | */ | |
862 | virtual TestF::Test* MakeTest() { return new Test(); } | |
863 | /** | |
864 | * Prefix on output | |
865 | * | |
866 | * @return The string "shift" | |
867 | */ | |
868 | virtual const char* Prefix() const { return "shift"; } | |
869 | /** | |
870 | * Make an output ntuple | |
871 | * | |
872 | * @return Possibly newly allocated NTuple | |
873 | */ | |
874 | virtual TNtuple* MakeNTuple () | |
875 | { | |
876 | return new TNtuple("shift","Shift scan", | |
877 | "Delta:xi:sigma:chi2nu:p:ep:c:ec"); | |
878 | } | |
879 | /** | |
880 | * Set-up for looping | |
881 | * | |
882 | * @param mode | |
883 | * @param nVal | |
884 | */ | |
885 | virtual void PreLoop(UShort_t mode, Int_t n) | |
886 | { | |
887 | if (mode == 0) { | |
888 | c2s = new TGraph2DErrors(n*n);c2s->SetName("cs"); | |
889 | p2s = new TGraph2DErrors(n*n);p2s->SetName("ps"); | |
890 | mcs = new TMultiGraph("mcs", "C as function of #sigma/#xi"); | |
891 | mps = new TMultiGraph("mps", "P as function of #sigma/#xi"); | |
892 | c2s->SetDirectory(0); | |
893 | p2s->SetDirectory(0); | |
894 | c2s->SetMarkerStyle(20); | |
895 | p2s->SetMarkerStyle(20); | |
896 | c2s->SetTitle("c(#xi,#sigma)"); | |
897 | p2s->SetTitle("p(#xi,#sigma)"); | |
898 | } | |
899 | else { | |
900 | const char* var = (mode == 1 ? "#sigma" : "#xi"); | |
901 | cs = new TGraphErrors(n); cs->SetName("cs"); | |
902 | ps = new TGraphErrors(n); ps->SetName("ps"); | |
903 | cs->SetMarkerStyle(20); | |
904 | ps->SetMarkerStyle(20); | |
905 | cs->SetTitle(Form("c(%s)", var)); | |
906 | ps->SetTitle(Form("p(%s)", var)); | |
907 | } | |
908 | Printf("%-8s | %-8s | %-8s | %-8s | %-8s %-9s | %-8s %-9s", | |
909 | "Delta_p", "xi", "sigma", "chi^2/nu", "c", "error", "p", "error"); | |
910 | } | |
911 | /** | |
912 | * Fill multi graph | |
913 | * | |
914 | * @param mg - multi-graph | |
915 | * @param i First index | |
916 | * @param j Second index | |
917 | * @param n Number of steps in each directon | |
918 | * @param u @f$\sigma/\xi@f$ | |
919 | * @param y Value | |
920 | * @param e Error | |
921 | */ | |
922 | void FillMG(TMultiGraph* mg, | |
923 | Int_t i, | |
924 | Int_t j, | |
925 | Int_t n, | |
926 | Double_t u, | |
927 | Double_t y, | |
928 | Double_t e) | |
929 | { | |
930 | TList* l = mg->GetListOfGraphs(); | |
931 | TGraphErrors* ge = (l ? static_cast<TGraphErrors*>(l->At(i)) : 0); | |
932 | if (!ge) { | |
933 | ge = new TGraphErrors(n); | |
934 | ge->SetName(Form("%s%02d", mg->GetName(), i)); | |
935 | ge->SetTitle(Form(Form("%d", i))); | |
936 | ge->SetMarkerColor(AliLandauGaus::GetIColor(i)); | |
937 | ge->SetLineColor(AliLandauGaus::GetIColor(i)); | |
938 | ge->SetMarkerStyle(20); | |
939 | mg->Add(ge); | |
940 | } | |
941 | ge->SetPoint(j, u, y); | |
942 | ge->SetPointError(j, 0, e); | |
943 | } | |
944 | /** | |
945 | * Process a single step | |
946 | * | |
947 | * @param mode Mode | |
948 | * @param i First index | |
949 | * @param j Second index | |
950 | * @param n Number of steps in each direction | |
951 | * @param bt Test result | |
952 | */ | |
953 | virtual void Step(UShort_t mode, | |
954 | Int_t i, | |
955 | Int_t j, | |
956 | Int_t n, | |
957 | TestF::Test* bt, | |
958 | TNtuple* nt) | |
959 | { | |
960 | Test* t = static_cast<Test*>(bt); | |
961 | Double_t xi = t->Xi(); | |
962 | Double_t sigma = t->Sigma(); | |
963 | Double_t c = t->fit->GetParameter(0); | |
964 | Double_t p = t->fit->GetParameter(1); | |
965 | Double_t ec = t->fit->GetParError(0); | |
966 | Double_t ep = t->fit->GetParError(1); | |
967 | Int_t idx = i * n + j; | |
968 | if (mode == 0) { | |
969 | c2s->SetPoint(idx, xi, sigma, c); | |
970 | p2s->SetPoint(idx, xi, sigma, p); | |
971 | c2s->SetPointError(idx, 0, 0, ec); | |
972 | p2s->SetPointError(idx, 0, 0, ep); | |
973 | FillMG(mcs, i, j, n, sigma/xi, c, ec); | |
974 | FillMG(mps, i, j, n, sigma/xi, p, ep); | |
975 | } | |
976 | else { | |
977 | cs->SetPoint(idx, (mode==1 ? sigma : xi), c); | |
978 | ps->SetPoint(idx, (mode==1 ? sigma : xi), p); | |
979 | cs->SetPointError(idx, 0, ec); | |
980 | ps->SetPointError(idx, 0, ep); | |
981 | } | |
982 | nt->Fill(1,xi,sigma,t->fit->GetChisquare()/t->fit->GetNDF(),c,ec,p,ep); | |
983 | } | |
984 | /** | |
985 | * Called after end of loop | |
986 | * | |
987 | * @param mode Mode of operation | |
988 | * @param can Canvas | |
989 | * @param nt NTuple | |
990 | * @param out Output file | |
991 | */ | |
992 | void PostLoop(UShort_t mode, TCanvas* can, TNtuple* /*nt*/, TFile* out) | |
993 | { | |
994 | can->Clear(); | |
995 | if (mode == 0) { | |
996 | can->Divide(2,2); | |
997 | can->cd(1); c2s->Draw("TRI2Z"); | |
998 | can->cd(2); p2s->Draw("TRI2Z"); | |
999 | can->cd(3); mcs->Draw("APL"); | |
1000 | can->cd(4); mps->Draw("APL"); | |
1001 | c2s->GetHistogram()->SetXTitle("#xi"); | |
1002 | p2s->GetHistogram()->SetXTitle("#xi"); | |
1003 | c2s->GetHistogram()->SetYTitle("#sigma"); | |
1004 | p2s->GetHistogram()->SetYTitle("#sigma"); | |
1005 | c2s->GetHistogram()->SetZTitle("c"); | |
1006 | p2s->GetHistogram()->SetZTitle("p"); | |
1007 | mcs->GetHistogram()->SetXTitle("#sigma/#xi"); | |
1008 | mps->GetHistogram()->SetXTitle("#sigma/#xi"); | |
1009 | mcs->GetHistogram()->SetYTitle("c"); | |
1010 | mps->GetHistogram()->SetYTitle("p"); | |
1011 | ||
1012 | out->cd(); | |
1013 | c2s->Write(); | |
1014 | p2s->Write(); | |
1015 | can->Write(); | |
1016 | mcs->Write(); | |
1017 | mps->Write(); | |
1018 | } | |
1019 | else { | |
1020 | can->Divide(2,1); | |
1021 | can->cd(1); cs->Draw("APL"); | |
1022 | can->cd(2); ps->Draw("APL"); | |
1023 | ||
1024 | out->cd(); | |
1025 | can->Write(); | |
1026 | cs->Write(); | |
1027 | ps->Write(); | |
1028 | } | |
1029 | PrintCanvas(can, 0, 0); | |
1030 | } | |
1031 | }; | |
1032 | #endif | |
1033 | ||
1034 | //==================================================================== | |
1035 | #ifdef TEST_FITTER | |
1036 | #include "AliLandauGausFitter.h" | |
1037 | #include <TRandom.h> | |
1038 | ||
1039 | struct TestFit : public TestF | |
1040 | { | |
1041 | struct Test : public TestF::Test | |
1042 | { | |
1043 | TH1* dist; | |
1044 | TF1* res; | |
1045 | TH1* pars; | |
1046 | Bool_t fDoShift; | |
1047 | Test() : TestF::Test(), dist(0), res(0), pars(0), fDoShift(true) {} | |
1048 | /** | |
1049 | * Re-implementation of TH1::FillRandom to accept a TF1 pointer | |
1050 | * | |
1051 | * @param h Histogram to fill | |
1052 | * @param f Function to use | |
1053 | * @param n Number of samples | |
1054 | */ | |
1055 | static void FillRandom(TH1* h, TF1* f, Int_t n=1000000) { | |
1056 | TAxis* xAxis = h->GetXaxis(); | |
1057 | ||
1058 | Int_t first = xAxis->GetFirst(); | |
1059 | Int_t last = xAxis->GetLast(); | |
1060 | Int_t nbinsx = last-first+1; | |
1061 | ||
1062 | TArrayD integral(nbinsx+1); | |
1063 | integral[0] = 0; | |
1064 | for (Int_t i = 1; i <= nbinsx; i++) { | |
1065 | Double_t fint = f->Integral(xAxis->GetBinLowEdge(i+first-1), | |
1066 | xAxis->GetBinUpEdge(i+first-1)); | |
1067 | integral[i] = integral[i-1] + fint; | |
1068 | } | |
1069 | ||
1070 | // - Normalize integral to 1 | |
1071 | if (integral[nbinsx] == 0) { | |
1072 | Error("FillRandom", "Integral = zero"); | |
1073 | return; | |
1074 | } | |
1075 | for (Int_t i = 1 ; i <= nbinsx; i++) | |
1076 | integral[i] /= integral[nbinsx]; | |
1077 | ||
1078 | // --------------Start main loop ntimes | |
1079 | for (Int_t i = 0; i < n; i++) { | |
1080 | Double_t r1 = gRandom->Rndm(i); | |
1081 | Int_t j = TMath::BinarySearch(nbinsx,&integral[0],r1); | |
1082 | Double_t x = (xAxis->GetBinLowEdge(j+first) + | |
1083 | xAxis->GetBinWidth(j+first) * | |
1084 | (r1-integral[j])/(integral[j+1] - integral[j])); | |
1085 | h->Fill(x); | |
1086 | } | |
1087 | } | |
1088 | /** | |
1089 | * Run the actual test | |
1090 | */ | |
1091 | void DoRun() | |
1092 | { | |
1093 | AliLandauGaus::EnableSigmaShift(1); // Enable sigma-shift | |
1094 | Double_t xMin = XMin(); | |
1095 | Double_t xMax = XMax(); | |
1096 | Int_t n = N(); | |
1097 | src->SetLineColor(kBlue); | |
1098 | src->SetLineStyle(2); | |
1099 | dist = new TH1D("dist", "Landau-Gaus distribution",500,xMin,xMax); | |
1100 | dist->SetXTitle("#Delta"); | |
1101 | dist->SetYTitle("dN/d#Delta"); | |
1102 | dist->SetFillStyle(3001); | |
1103 | dist->SetFillColor(kMagenta+1); | |
1104 | dist->Sumw2(); | |
1105 | dist->SetDirectory(0); | |
1106 | FillRandom(dist, src); | |
1107 | ||
1108 | Int_t peakBin = dist->GetMaximumBin(); | |
1109 | Double_t max = dist->GetBinContent(peakBin); | |
1110 | dist->Scale(1/max); | |
1111 | ||
1112 | AliLandauGaus::EnableSigmaShift(fDoShift); // Disable sigma-shift | |
1113 | AliLandauGausFitter f(xMin, xMax, 10); | |
1114 | // f.SetDebug(true); | |
1115 | ||
1116 | TF1* r = f.FitNParticle(dist, n); | |
1117 | if (!r) return; | |
1118 | r->SetName("fit"); | |
1119 | ||
1120 | // Get the status for the covariance calculation: | |
1121 | // 0: Not calculatuted | |
1122 | // 1: Approximate | |
1123 | // 2: Forced positive definite | |
1124 | // 3: Accurate | |
1125 | Int_t covStatus = | |
1126 | static_cast<TFitResult*>(f.GetFitResults().At(n-1))->CovMatrixStatus(); | |
1127 | if (covStatus < 2) ::Warning("", "Fit may not be good"); | |
1128 | ||
1129 | r->SetRange(xMin, xMax); | |
1130 | res = static_cast<TF1*>(r->Clone()); | |
1131 | res->SetLineColor(kPink+10); | |
1132 | res->SetLineStyle(1); | |
1133 | res->SetLineWidth(3); | |
1134 | dist->GetListOfFunctions()->Add(res); | |
1135 | ||
1136 | Double_t min = 0.1*dist->GetMinimum(); | |
1137 | TF1* old = src; // Tuck away src | |
1138 | src = res; // Use fit now | |
1139 | for (Int_t i = 1; i <= n; i++) { | |
1140 | Double_t rDelta = Delta(); | |
1141 | Double_t rxi = Xi(); | |
1142 | Double_t rsigma = Sigma(); | |
1143 | TF1* comp = Fi(i); | |
1144 | comp->SetLineWidth(2); | |
1145 | comp->SetLineStyle(1); | |
1146 | ||
1147 | // Modifies arguments! | |
1148 | AliLandauGaus::IPars(i, rDelta, rxi, rsigma); | |
1149 | TLine* l = VerticalLine(1, comp->GetLineColor(), rDelta, min, | |
1150 | 1.1*old->Eval(rDelta)); | |
1151 | dist->GetListOfFunctions()->Add(comp); | |
1152 | dist->GetListOfFunctions()->Add(l); | |
1153 | } | |
1154 | src = old; // restore | |
1155 | ||
1156 | } | |
1157 | /** | |
1158 | * Draw results in pad | |
1159 | * | |
1160 | * @param p Pad | |
1161 | */ | |
1162 | void DrawInPad(TVirtualPad* p) | |
1163 | { | |
1164 | // Double_t scale = src->GetMaximum(); | |
1165 | // src->SetParameter(0, 1/scale); | |
1166 | p->Clear(); | |
1167 | p->SetTopMargin(0.0); | |
1168 | p->SetRightMargin(0.0); | |
1169 | p->SetBottomMargin(0.0); | |
1170 | p->SetLeftMargin(0.0); | |
1171 | p->SetFillColor(kGray); | |
1172 | p->Divide(2,1,0.001,0.001); | |
1173 | TVirtualPad* q = p->cd(1); | |
1174 | q->SetLogy(); | |
1175 | q->SetTopMargin(0.01); | |
1176 | q->SetRightMargin(0.02); | |
1177 | q->SetFillColor(kWhite); | |
1178 | gStyle->SetOptStat(0); | |
1179 | gStyle->SetOptFit(11111); | |
1180 | gStyle->SetOptTitle(false); | |
1181 | gStyle->SetStatY(1-q->GetTopMargin()); | |
1182 | ||
1183 | dist->DrawCopy("hist"); | |
1184 | dist->DrawCopy("e same"); | |
1185 | AliLandauGaus::EnableSigmaShift(1); // Enable sigma-shift | |
1186 | src->DrawClone("same"); | |
1187 | AliLandauGaus::EnableSigmaShift(fDoShift); // Disable sigma-shift | |
1188 | if (!res) return; | |
1189 | ||
1190 | ||
1191 | q = p->cd(2); | |
1192 | q->SetTopMargin(0.01); | |
1193 | q->SetRightMargin(0.02); | |
1194 | q->SetFillColor(kWhite); | |
1195 | Int_t nPar = src->GetNpar(); | |
1196 | pars = new TH1D("pars", "#chi^{2}/#nu & #Deltap_{i}/#deltap_{i}", | |
1197 | nPar+1, 0, nPar+1); | |
1198 | pars->SetDirectory(0); | |
1199 | pars->SetFillColor(kCyan+2); | |
1200 | pars->SetFillStyle(3001); | |
1201 | pars->SetYTitle("#Deltap_{i}/#deltap_{i}"); | |
1202 | ||
1203 | TPaveStats* stats = new TPaveStats(0,0,0,0); | |
1204 | stats->SetX1NDC(gStyle->GetStatX()-gStyle->GetStatW()); | |
1205 | stats->SetY1NDC(gStyle->GetStatY()-gStyle->GetStatH()); | |
1206 | stats->SetX2NDC(gStyle->GetStatX()); | |
1207 | stats->SetY2NDC(gStyle->GetStatY()); | |
1208 | stats->SetTextFont(42); | |
1209 | ||
1210 | stats->SetBorderSize(1); | |
1211 | stats->SetFillColor(kWhite); | |
1212 | stats->AddText(""); | |
1213 | pars->GetListOfFunctions()->Add(stats); | |
1214 | TString fmt(Form("%%s = %%%s", stats->GetFitFormat())); | |
1215 | for (Int_t i = 0; i < nPar; i++) { | |
1216 | Double_t low, high; | |
1217 | res->GetParLimits(i, low, high); | |
1218 | Double_t pSrc = src->GetParameter(i); | |
1219 | Double_t pRes = res->GetParameter(i); | |
1220 | Double_t eRes = res->GetParError(i); | |
1221 | Bool_t fix = (low == high && ((pRes == 0 && low == 1) || low==pRes)); | |
1222 | Double_t diff = fix ? 0 : TMath::Abs(pSrc-pRes)/eRes; | |
1223 | pars->GetXaxis()->SetBinLabel(i+1,res->GetParName(i)); | |
1224 | pars->SetBinContent(i+1, diff); | |
1225 | stats->AddText(Form(fmt.Data(), res->GetParName(i), pSrc)); | |
1226 | } | |
1227 | pars->GetXaxis()->SetBinLabel(nPar+1,"#chi^{2}/#nu"); | |
1228 | pars->SetBinContent(nPar+1,res->GetChisquare() / res->GetNDF()); | |
1229 | pars->DrawCopy(); | |
1230 | q->Modified(); | |
1231 | q->Update(); | |
1232 | ||
1233 | p->Modified(); | |
1234 | p->Update(); | |
1235 | p->cd(); | |
1236 | ||
1237 | PrintFs(false); | |
1238 | } | |
1239 | /** | |
1240 | * Print results | |
1241 | * | |
1242 | * @param full | |
1243 | */ | |
1244 | void PrintFs(Bool_t full=false) { | |
1245 | if (full) { | |
1246 | Int_t nPar = src->GetNpar(); | |
1247 | Printf("%-2s %-10s | %-8s | %8s %-9s | %9.5s", | |
1248 | "#", "Name", "Source", "Fit", "Error", "delta"); | |
1249 | for (Int_t i = 0; i < nPar; i++) { | |
1250 | Double_t low, high; | |
1251 | res->GetParLimits(i, low, high); | |
1252 | Double_t pSrc = src->GetParameter(i); | |
1253 | Double_t pRes = res->GetParameter(i); | |
1254 | Double_t eRes = res->GetParError(i); | |
1255 | Bool_t fix = (low == high && ((pRes == 0 && low==1) || low==pRes)); | |
1256 | Double_t diff = fix ? 0 : TMath::Abs(pSrc-pRes)/eRes; | |
1257 | Printf("%2d %10s | %8.4f | %8.4f %9.5f | %9.5f %s", | |
1258 | i, src->GetParName(i), pSrc, pRes, eRes, diff, | |
1259 | (fix ? "fixed" : (diff < 2 ? "ok" : "bad"))); | |
1260 | } | |
1261 | Printf("chi^2/nu = %8.4f/%-3d = %9.5f", | |
1262 | res->GetChisquare(), res->GetNDF(), | |
1263 | res->GetChisquare() / res->GetNDF()); | |
1264 | } | |
1265 | else | |
1266 | Printf("%8.4f", res->GetChisquare() / res->GetNDF()); | |
1267 | } | |
1268 | /** | |
1269 | * Write results to disk | |
1270 | * | |
1271 | * @param d Directory to write to | |
1272 | */ | |
1273 | virtual void WriteOut(TDirectory* d) | |
1274 | { | |
1275 | d->cd(); | |
1276 | src->Write(); | |
1277 | dist->Write(); | |
1278 | if (res) res->Write(); | |
1279 | if (pars) pars->Write(); | |
1280 | } | |
1281 | ||
1282 | }; | |
1283 | //__________________________________________________________________ | |
1284 | TestF::Test* MakeTest() { return new Test(); } | |
1285 | const char* Prefix() const { return "fit"; } | |
1286 | void PreLoop(UShort_t,Int_t) { | |
1287 | Printf("%-8s | %-8s | %-8s | %-8s ", "Delta_p", "xi", "sigma", "chi^2/nu"); | |
1288 | } | |
1289 | void Step(UShort_t,Int_t,Int_t,Int_t,TestF::Test*,TNtuple*) {} | |
1290 | void PostLoop(UShort_t,TCanvas* c,TNtuple*,TFile*) { PrintCanvas(c,0,0); } | |
1291 | }; | |
1292 | #endif | |
1293 | // | |
1294 | // EOF | |
1295 | // | |
1296 | ||
1297 |