AliPhysics  251aa1e (251aa1e)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TestF.C
Go to the documentation of this file.
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 //====================================================================
25 struct TestF
26 {
27  //__________________________________________________________________
31  struct Test {
33  TF1* src;
37  Test() : src(0) {}
41  Test(const Test& o)
42  : src(static_cast<TF1*>(o.src ? o.src->Clone() : 0)) {}
46  virtual ~Test() { delete src; }
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  }
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  }
94  Double_t C() const { return src->GetParameter(AliLandauGaus::kC); }
98  Double_t Delta() const { return src->GetParameter(AliLandauGaus::kDelta); }
102  Double_t Xi() const { return src->GetParameter(AliLandauGaus::kXi); }
106  Double_t Sigma() const { return src->GetParameter(AliLandauGaus::kSigma); }
110  Double_t SigmaN() const { return src->GetParameter(AliLandauGaus::kSigmaN);}
114  Double_t N() const { return src->GetParameter(AliLandauGaus::kN); }
118  Double_t A(Int_t i) const {
119  return (i <= 1 ? 1 : src->GetParameter(AliLandauGaus::kN+(i-1))); }
123  Double_t XMin() const { return src->GetXmin(); }
127  Double_t XMax() const { return src->GetXmax(); }
141  Double_t sigmaN, Int_t i) const {
142  return AliLandauGaus::MakeFi(c,delta,xi,sigma,sigmaN,i,XMin(),XMax());
143  }
152  TF1* Fi(Int_t i) {
153  return Fi(C()*A(i),Delta(),Xi(),Sigma(),SigmaN(),i);
154  }
155 
163  virtual void DoRun() = 0;
169  virtual void WriteOut(TDirectory* d) = 0;
175  virtual void DrawInPad(TVirtualPad* p) = 0;
176  /* @} */
177  };
178  //__________________________________________________________________
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() {}
209  virtual Test* MakeTest() = 0;
225  Test* DoOne(const Double_t c,
226  const Double_t delta,
227  const Double_t xi,
228  const Double_t sigma,
229  const Int_t n,
230  const Double_t* a,
231  const Double_t xMin,
232  const Double_t xMax,
233  TVirtualPad* pad=0)
234  {
235  Test* t = MakeTest();
236  t->Run(c, delta, xi, sigma, n, a, xMin, xMax);
237 
238  TString name(Form("sigma%04d_xi%04d", Int_t(100*sigma), Int_t(100*xi)));
239  TString title(Form("Sigma=%f xi=%f", sigma, xi));
240  if (!pad) {
241  TCanvas* can = new TCanvas(name, title, 1200, 900);
242  pad = can;
243  }
244  else
245  pad->Clear();
246  t->DrawInPad(pad);
247 
248  return t;
249  }
261  const Double_t xi=1,
262  const Int_t n=10,
263  TVirtualPad* pad=0)
264  {
265  const Double_t c = 1;
266  const Double_t delta = 10;
267  Double_t idelta = delta;
268  Double_t ixi = xi;
269  Double_t isigma = sigma;
270  AliLandauGaus::IPars(n+1,idelta,ixi, isigma);
271  Double_t up = idelta + ixi + isigma; // 8*n
272  Double_t a[n];
273  a[0] = 1;
274  for (Int_t i = 1; i < n; i++)
275  a[i] = TMath::Exp(0.5 - TMath::Sqrt(2)*(i+1));
276 
277  return DoOne(c,delta,xi,sigma,n,a,-3,up, pad);
278  }
286  Test* DoRealistic(TVirtualPad* pad=0)
287  {
288  Double_t c = 0.4070;
289  Double_t delta = 0.5651;
290  Double_t xi = 0.0643;
291  Double_t sigma = 0.0611;
292  Int_t n = 5;
293  Double_t a[] = {1, 0.1179, 0.0333, 0.0068, 0.0012 };
294 
295  return DoOne(c, delta, xi, sigma, n, a, 0.02, n, pad);
296  }
303  virtual const char* Prefix() const = 0;
309  virtual TNtuple* MakeNTuple() { return 0; }
318  virtual void SetupScan(UShort_t mode,
319  TFile*& file,
320  TCanvas*& canvas,
321  TNtuple*& nt)
322  {
323  TString base = Form("%s%s", Prefix(),
324  (mode == 0 ? "SigmaXi" :
325  mode == 1 ? "Sigma" : "Xi"));
326  file = TFile::Open(Form("%s.root", base.Data()), "RECREATE");
327  nt = MakeNTuple();
328  canvas = new TCanvas(base, Form("%s.pdf", base.Data()),
329  TMath::Sqrt(2)*900, 900);
330  TString title("pdf Landscape Title:Scan");
331  gSystem->RedirectOutput("/dev/null");
332  canvas->Print(Form("%s[", canvas->GetTitle()), title);
333  gSystem->RedirectOutput(0);
334  }
344  void PrintCanvas(TCanvas* c, Double_t xi, Double_t sigma) {
345  TString tit = "pdf Landscape Title:";
346  Bool_t cls = false;
347  if (xi > 0 && sigma > 0) {
348  tit.Append(Form("xi=%8.4f sigma=%8.4f", xi, sigma));
349  }
350  else {
351  tit.Append("Results");
352  cls = true;
353  }
354  gSystem->RedirectOutput("/dev/null");
355  c->Print(c->GetTitle(), tit);
356  if (cls) {
357  c->Clear();
358  c->Print(Form("%s]", c->GetTitle()), tit);
359  }
360  gSystem->RedirectOutput(0);
361  if (cls)
362  ::Info("", "Plots saved in %s", c->GetTitle());
363  }
370  virtual void PreLoop(UShort_t mode, Int_t nVal) = 0;
381  virtual void Step(UShort_t mode,
382  Int_t i,
383  Int_t j,
384  Int_t n,
385  Test* test,
386  TNtuple* nt) = 0;
395  virtual void PostLoop(UShort_t mode,
396  TCanvas* can,
397  TNtuple* nt,
398  TFile* out) = 0;
407  void ScanOne(Bool_t scanSigma,
408  Int_t n,
409  const Double_t* values,
410  Int_t maxN=10)
411  {
412  UShort_t mod = scanSigma ? 1 : 2;
413  Scan(mod, n, values, maxN);
414  }
422  void ScanOne(Bool_t scanSigma,
423  const TArrayD& values,
424  Int_t maxN)
425  {
426  ScanOne(scanSigma, values.GetSize(), values.GetArray(), maxN);
427  }
435  void ScanTwo(Int_t n,
436  const Double_t* values,
437  Int_t maxN=10)
438  {
439  Scan(0, n, values, maxN);
440  }
447  void ScanTwo(const TArrayD& values, Int_t maxN)
448  {
449  ScanTwo(values.GetSize(), values.GetArray(), maxN);
450  }
460  Int_t n,
461  const Double_t* values,
462  Int_t maxN)
463  {
464  TFile* out = 0;
465  TNtuple* nt = 0;
466  TCanvas* can = 0;
467  SetupScan(mode, out, can, nt);
468  PreLoop(mode, n);
469 
470  Test* t = DoRealistic(can);
471  t->WriteOut(out->mkdir("realistic"));
472  can->Write();
473  PrintCanvas(can, t->Xi(), t->Sigma());
474 
475  if (mode == 0) LoopTwo(n, values, maxN, can, nt, out);
476  else LoopOne(mode, n, values, maxN, can, nt, out);
477 
478  PostLoop(mode, can, nt, out);
479  out->Write();
480  out->Close();
481 
482  }
495  Int_t n,
496  const Double_t* values,
497  Int_t maxN,
498  TCanvas* can,
499  TNtuple* nt,
500  TFile* out) {
501  Bool_t scanSigma = mode == 1;
502  const char* var = (scanSigma ? "sigma" : "xi");
503 
504  for (Int_t i = 0; i < n; i++) {
505  Double_t v = values[i];
506  Double_t sigma = (scanSigma ? v : 1);
507  Double_t xi = (!scanSigma ? v : 1);
508 
509  StepIt(mode, xi, sigma, maxN, 0, i, n, can, nt, out, var, v);
510  }
511  }
522  void LoopTwo(Int_t n,
523  const Double_t* values,
524  Int_t maxN,
525  TCanvas* can,
526  TNtuple* nt,
527  TFile* out)
528  {
529  UShort_t mode = 0;
530  for (Int_t i = 0; i < n; i++) {
531  Double_t xi = values[i];
532  TDirectory* d = out->mkdir(Form("xi%04d", Int_t(100*xi)));
533  for (Int_t j = 0; j < n; j++) {
534  Double_t sigma = values[j];
535  StepIt(mode, xi, sigma, maxN, i, j, n, can, nt, d, "sigma", sigma);
536  }
537  }
538  }
556  Double_t xi,
557  Double_t sigma,
558  Int_t maxN,
559  Int_t i,
560  Int_t j,
561  Int_t n,
562  TCanvas* can,
563  TNtuple* nt,
564  TDirectory* p,
565  const char* pre,
566  Double_t v)
567  {
568  TDirectory* d = p->mkdir(Form("%s%04d", pre, Int_t(100*v)));
569  Test* t = DoUnit(sigma, xi, maxN, can);
570  t->WriteOut(d);
571  can->Write();
572  PrintCanvas(can, xi, sigma);
573 
574  Step(mode, i, j, n, t, nt);
575  }
576 };
577 
578 //====================================================================
579 #ifdef TEST_SHIFT
580 #include "WrappedGraph.C"
581 #include <TGraphErrors.h>
582 #include <TGraph2DErrors.h>
583 #include <TMultiGraph.h>
584 #include <Math/RootFinder.h>
585 
586 struct TestShift : public TestF
587 {
604  static Double_t fitFunc(Double_t* xp, Double_t* pp) {
605  Double_t x = *xp;
606  Double_t c = pp[0];
607  Double_t p = pp[1];
608  Double_t xi = pp[2];
609  Double_t sigma = pp[3];
610  Double_t u = sigma / xi;
611  Double_t a = c * u * sigma;
612  Double_t b = p * u * TMath::Sqrt(u);
613 
614  return a / TMath::Power(1+1./x, b);
615  }
622  static void ScaleGraph(TGraph* g, Double_t scale)
623  {
624  // Info("", "Scaling graph %p with %f", g, scale);
625  if (scale == 1) return;
626  for (Int_t i = 0; i < g->GetN(); i++)
627  g->SetPoint(i,g->GetX()[i], g->GetY()[i] * scale);
628  }
629  //__________________________________________________________________
630  struct Test : public TestF::Test
631  {
632  TGraph* zeros; // Zeros of the derivatives
633  TGraph* diff; // Difference from calculated Delta to roots
634  TList comps; // Components
635  TList derivs; // Derivatives
636  TF1* fit; // Fit to difference
637  Test() : TestF::Test(), zeros(0), diff(0), fit(0)
638  {
639  comps.SetName("components");
640  derivs.SetName("derivatives");
641 
642  }
646  void DoRun()
647  {
648  AliLandauGaus::EnableSigmaShift(false); // Disable sigma-shift
649  const Int_t n = N();
650  zeros = new TGraph(n);
651  zeros->SetName("zeros");
652  zeros->SetTitle("Zeros of derivatives");
653  zeros->SetMarkerStyle(25);
654  zeros->SetMarkerSize(1.2);
655  zeros->SetMarkerColor(kRed+1);
656 
657  diff = new TGraph(n);
658  diff->SetName("diff");
659  diff->SetTitle("#delta#Delta_{p}");
660  diff->SetMarkerStyle(20);
661  diff->SetMarkerColor(kBlue+1);
662 
663  // Make the components
664  Int_t j = 0;
665  for (Int_t i = 1; i <= n; i++, j++)
666  if (!MakeComp(i)) break;
667  zeros->Set(j);
668  diff->Set(j);
669 
670  // Now fit
671  fit = new TF1("fit", &fitFunc, 1, j, 4);
672  fit->SetParNames("c", "p", "#xi", "#sigma");
673  fit->SetParameter(0, .5);
674  fit->SetParameter(1, .5);
675  fit->FixParameter(2, Xi());
676  fit->FixParameter(3, Sigma());
677 
678  // fit->SetParNames("c", "#xi", "#sigma");
679  // fit->FixParameter(1, xi);
680 
681 
682  Int_t ifit = 4;
683  do {
684  gSystem->RedirectOutput("/dev/null");
685  diff->Fit(fit, Form("MR%s", ifit == 0 ? "Q" : "Q"));
686  gSystem->RedirectOutput(0);
687  ifit--;
688  } while (ifit >= 0);
689  Printf("%8.4f | %8.4f %9.5f | %8.4f %9.5f",
690  fit->GetChisquare() / fit->GetNDF(),
691  fit->GetParameter(0), fit->GetParError(0),
692  fit->GetParameter(1), fit->GetParError(1));
693 
694  }
700  Bool_t MakeComp(const Int_t i) {
701  TF1* fi = Fi(i);
702  fi->SetNpx(src->GetNpx());
703  // fi->Print();
704 
705  TGraph* dF = new TGraph(fi, "d");
706  dF->SetName(Form("d%s", fi->GetName()));
707  dF->SetTitle(Form("d%s/d#Delta", fi->GetTitle()));
708  dF->SetMarkerStyle(1); // A dot
709  // dF->SetMarkerStyle(20+i-1);
710  dF->SetMarkerColor(fi->GetLineColor());
711  dF->SetLineColor(fi->GetLineColor());
712 
713  Double_t max = TMath::MaxElement(dF->GetN(), dF->GetY());
714  if (max < 1e-6) return false;
715 
716  ScaleGraph(dF, 1./max);
717  Double_t delta = Delta();
718  Double_t xi = Xi();
719  Double_t deltaI = i * (delta + xi * TMath::Log(i));
720  Double_t xiI = i * xi;
721  Double_t xR = FindRoot(dF, deltaI-2*xiI, deltaI+2*xiI);
722 
723  // Printf("Component %2d: Exp: %8.4f Got: %8.4f Diff: %9.5f",
724  // i, deltaI, xR, (xR-deltaI));
725  // printf(".");
726  comps.Add(fi);
727  derivs.Add(dF);
728  zeros->SetPoint(i-1, xR, dF->Eval(xR));
729  diff->SetPoint(i-1, i, xR - deltaI);
730 
731  return true;
732  }
742  Double_t FindRoot(TGraph* g, Double_t xMin, Double_t xMax) {
743  WrappedGraph* wG = new WrappedGraph(g);
744  ROOT::Math::RootFinder rfb(ROOT::Math::RootFinder::kBRENT);
745  rfb.SetFunction(*wG, xMin, xMax);
746  rfb.Solve();
747  Double_t xR = rfb.Root();
748  return xR;
749  }
755  void DrawInPad(TVirtualPad* body) {
756  body->SetTopMargin(0.01);
757  body->SetRightMargin(0.01);
758  body->Divide(2,1);
759 
760  TVirtualPad* p = body->cd(1);
761  p->SetTopMargin(0.01);
762  p->SetRightMargin(0.01);
763  p->Divide(1,2,0,0);
764 
765  TVirtualPad* q = p->cd(1);
766  q->SetRightMargin(0.01);
767  q->SetGridx();
768  q->SetLogy();
769  src->Draw();
770  Double_t fmax = src->GetHistogram()->GetMaximum();
771  Double_t fmin = src->GetHistogram()->GetMinimum();
772  fmin = TMath::Max(fmin, fmax/1e8);
773  src->GetHistogram()->SetMinimum(fmin);
774  TLatex* ltx = new TLatex(q->GetLeftMargin()+.02,
775  1-q->GetTopMargin()-.02,
776  Form("#Delta_{p}=%f, #xi=%f, #sigma=%f",
777  Delta(), Xi(), Sigma()));
778  // Printf("%s: Max=%f scale=%f", src->GetName(), fmax, 1/fmax);
779  ltx->SetTextFont(42);
780  ltx->SetNDC();
781  ltx->SetTextAlign(13);
782  ltx->Draw();
783 
784  q = p->cd(2);
785  q->SetRightMargin(0.01);
786  q->SetGridx();
787 
788  Double_t x = 1-q->GetRightMargin()-.02;
789  Double_t y = 1-q->GetTopMargin()-.1;
790  for (Int_t i = 1; i <= comps.GetEntries(); i++) {
791  p->cd(1);
792  TF1* fi = static_cast<TF1*>(comps.At(i-1));
793  fi->Draw("same");
794 
795  Double_t deltaR = zeros->GetX()[i-1];
796  Double_t deltaI = i * (Delta() + Xi() * TMath::Log(i));
797 
798  Color_t col = fi->GetLineColor();
799  VerticalLine(1,col,deltaI,fmin,fmax)->Draw();
800  VerticalLine(3,col,deltaR,fmin,fmax)->Draw();
801 
802  ltx = new TLatex(x,y,Form("#Delta_{p,%2d}=%8.4f %8.4f",
803  i,deltaI,deltaR));
804  ltx->SetTextAlign(33);
805  ltx->SetTextFont(42);
806  ltx->SetNDC();
807  ltx->Draw();
808  y -= ltx->GetTextSize() + .01;
809 
810  p->cd(2);
811 
812  TGraph* df = static_cast<TGraph*>(derivs.At(i-1));
813  df->Draw(i == 1 ? "alp" : "lp");
814  df->GetHistogram()->GetXaxis()->SetRangeUser(XMin(), XMax());
815  VerticalLine(1,col,deltaI,-.6,1.1)->Draw();
816  VerticalLine(3,col,deltaR,-.6,1.1)->Draw();
817  }
818  p->cd(2);
819  zeros->Draw("p");
820 
821  gStyle->SetOptFit(111111);
822  gStyle->SetStatY(.6);
823  p = body->cd(2);
824  p->SetRightMargin(0.02);
825  p->SetTopMargin(0.02);
826  diff->Draw("alp");
827  diff->GetHistogram()->GetListOfFunctions()->Add(fit);
828  diff->GetHistogram()->SetXTitle("N_{particle}");
829  diff->GetHistogram()->SetYTitle("#Delta_{r}-#Delta_{p}");
830  p->Clear();
831  diff->Draw("alp");
832 
833  body->Modified();
834  body->Update();
835  body->cd();
836  }
842  void WriteOut(TDirectory* dir) {
843  dir->cd();
844  src->Write();
845  zeros->Write();
846  diff->Write();
847  comps.Write(comps.GetName(), TObject::kSingleKey);
848  derivs.Write(derivs.GetName(), TObject::kSingleKey);
849  fit->Write();
850  }
851  };
852  //__________________________________________________________________
853  TGraphErrors* cs;
854  TGraphErrors* ps;
855  TGraph2DErrors* c2s;
856  TGraph2DErrors* p2s;
857  TMultiGraph* mcs;
858  TMultiGraph* mps;
859 
865  virtual TestF::Test* MakeTest() { return new Test(); }
871  virtual const char* Prefix() const { return "shift"; }
877  virtual TNtuple* MakeNTuple ()
878  {
879  return new TNtuple("shift","Shift scan",
880  "Delta:xi:sigma:chi2nu:p:ep:c:ec");
881  }
888  virtual void PreLoop(UShort_t mode, Int_t n)
889  {
890  if (mode == 0) {
891  c2s = new TGraph2DErrors(n*n);c2s->SetName("cs");
892  p2s = new TGraph2DErrors(n*n);p2s->SetName("ps");
893  mcs = new TMultiGraph("mcs", "C as function of #sigma/#xi");
894  mps = new TMultiGraph("mps", "P as function of #sigma/#xi");
895  c2s->SetDirectory(0);
896  p2s->SetDirectory(0);
897  c2s->SetMarkerStyle(20);
898  p2s->SetMarkerStyle(20);
899  c2s->SetTitle("c(#xi,#sigma)");
900  p2s->SetTitle("p(#xi,#sigma)");
901  }
902  else {
903  const char* var = (mode == 1 ? "#sigma" : "#xi");
904  cs = new TGraphErrors(n); cs->SetName("cs");
905  ps = new TGraphErrors(n); ps->SetName("ps");
906  cs->SetMarkerStyle(20);
907  ps->SetMarkerStyle(20);
908  cs->SetTitle(Form("c(%s)", var));
909  ps->SetTitle(Form("p(%s)", var));
910  }
911  Printf("%-8s | %-8s | %-8s | %-8s | %-8s %-9s | %-8s %-9s",
912  "Delta_p", "xi", "sigma", "chi^2/nu", "c", "error", "p", "error");
913  }
925  void FillMG(TMultiGraph* mg,
926  Int_t i,
927  Int_t j,
928  Int_t n,
929  Double_t u,
930  Double_t y,
931  Double_t e)
932  {
933  TList* l = mg->GetListOfGraphs();
934  TGraphErrors* ge = (l ? static_cast<TGraphErrors*>(l->At(i)) : 0);
935  if (!ge) {
936  ge = new TGraphErrors(n);
937  ge->SetName(Form("%s%02d", mg->GetName(), i));
938  ge->SetTitle(Form(Form("%d", i)));
939  ge->SetMarkerColor(AliLandauGaus::GetIColor(i));
940  ge->SetLineColor(AliLandauGaus::GetIColor(i));
941  ge->SetMarkerStyle(20);
942  mg->Add(ge);
943  }
944  ge->SetPoint(j, u, y);
945  ge->SetPointError(j, 0, e);
946  }
956  virtual void Step(UShort_t mode,
957  Int_t i,
958  Int_t j,
959  Int_t n,
960  TestF::Test* bt,
961  TNtuple* nt)
962  {
963  Test* t = static_cast<Test*>(bt);
964  Double_t xi = t->Xi();
965  Double_t sigma = t->Sigma();
966  Double_t c = t->fit->GetParameter(0);
967  Double_t p = t->fit->GetParameter(1);
968  Double_t ec = t->fit->GetParError(0);
969  Double_t ep = t->fit->GetParError(1);
970  Int_t idx = i * n + j;
971  if (mode == 0) {
972  c2s->SetPoint(idx, xi, sigma, c);
973  p2s->SetPoint(idx, xi, sigma, p);
974  c2s->SetPointError(idx, 0, 0, ec);
975  p2s->SetPointError(idx, 0, 0, ep);
976  FillMG(mcs, i, j, n, sigma/xi, c, ec);
977  FillMG(mps, i, j, n, sigma/xi, p, ep);
978  }
979  else {
980  cs->SetPoint(idx, (mode==1 ? sigma : xi), c);
981  ps->SetPoint(idx, (mode==1 ? sigma : xi), p);
982  cs->SetPointError(idx, 0, ec);
983  ps->SetPointError(idx, 0, ep);
984  }
985  nt->Fill(1,xi,sigma,t->fit->GetChisquare()/t->fit->GetNDF(),c,ec,p,ep);
986  }
995  void PostLoop(UShort_t mode, TCanvas* can, TNtuple* /*nt*/, TFile* out)
996  {
997  can->Clear();
998  if (mode == 0) {
999  can->Divide(2,2);
1000  can->cd(1); c2s->Draw("TRI2Z");
1001  can->cd(2); p2s->Draw("TRI2Z");
1002  can->cd(3); mcs->Draw("APL");
1003  can->cd(4); mps->Draw("APL");
1004  c2s->GetHistogram()->SetXTitle("#xi");
1005  p2s->GetHistogram()->SetXTitle("#xi");
1006  c2s->GetHistogram()->SetYTitle("#sigma");
1007  p2s->GetHistogram()->SetYTitle("#sigma");
1008  c2s->GetHistogram()->SetZTitle("c");
1009  p2s->GetHistogram()->SetZTitle("p");
1010  mcs->GetHistogram()->SetXTitle("#sigma/#xi");
1011  mps->GetHistogram()->SetXTitle("#sigma/#xi");
1012  mcs->GetHistogram()->SetYTitle("c");
1013  mps->GetHistogram()->SetYTitle("p");
1014 
1015  out->cd();
1016  c2s->Write();
1017  p2s->Write();
1018  can->Write();
1019  mcs->Write();
1020  mps->Write();
1021  }
1022  else {
1023  can->Divide(2,1);
1024  can->cd(1); cs->Draw("APL");
1025  can->cd(2); ps->Draw("APL");
1026 
1027  out->cd();
1028  can->Write();
1029  cs->Write();
1030  ps->Write();
1031  }
1032  PrintCanvas(can, 0, 0);
1033  }
1034 };
1035 #endif
1036 
1037 //====================================================================
1038 #ifdef TEST_FITTER
1039 #include "AliLandauGausFitter.h"
1040 #include <TRandom.h>
1041 
1042 struct TestFit : public TestF
1043 {
1044  struct Test : public TestF::Test
1045  {
1046  TH1* dist;
1047  TF1* res;
1048  TH1* pars;
1049  Bool_t fDoShift;
1050  Test() : TestF::Test(), dist(0), res(0), pars(0), fDoShift(true) {}
1058  static void FillRandom(TH1* h, TF1* f, Int_t n=1000000) {
1059  TAxis* xAxis = h->GetXaxis();
1060 
1061  Int_t first = xAxis->GetFirst();
1062  Int_t last = xAxis->GetLast();
1063  Int_t nbinsx = last-first+1;
1064 
1065  TArrayD integral(nbinsx+1);
1066  integral[0] = 0;
1067  for (Int_t i = 1; i <= nbinsx; i++) {
1068  Double_t fint = f->Integral(xAxis->GetBinLowEdge(i+first-1),
1069  xAxis->GetBinUpEdge(i+first-1));
1070  integral[i] = integral[i-1] + fint;
1071  }
1072 
1073  // - Normalize integral to 1
1074  if (integral[nbinsx] == 0) {
1075  Error("FillRandom", "Integral = zero");
1076  return;
1077  }
1078  for (Int_t i = 1 ; i <= nbinsx; i++)
1079  integral[i] /= integral[nbinsx];
1080 
1081  // --------------Start main loop ntimes
1082  for (Int_t i = 0; i < n; i++) {
1083  Double_t r1 = gRandom->Rndm(i);
1084  Int_t j = TMath::BinarySearch(nbinsx,&integral[0],r1);
1085  Double_t x = (xAxis->GetBinLowEdge(j+first) +
1086  xAxis->GetBinWidth(j+first) *
1087  (r1-integral[j])/(integral[j+1] - integral[j]));
1088  h->Fill(x);
1089  }
1090  }
1094  void DoRun()
1095  {
1096  AliLandauGaus::EnableSigmaShift(1); // Enable sigma-shift
1097  Double_t xMin = XMin();
1098  Double_t xMax = XMax();
1099  Int_t n = N();
1100  src->SetLineColor(kBlue);
1101  src->SetLineStyle(2);
1102  dist = new TH1D("dist", "Landau-Gaus distribution",500,xMin,xMax);
1103  dist->SetXTitle("#Delta");
1104  dist->SetYTitle("dN/d#Delta");
1105  dist->SetFillStyle(3001);
1106  dist->SetFillColor(kMagenta+1);
1107  dist->Sumw2();
1108  dist->SetDirectory(0);
1109  FillRandom(dist, src);
1110 
1111  Int_t peakBin = dist->GetMaximumBin();
1112  Double_t max = dist->GetBinContent(peakBin);
1113  dist->Scale(1/max);
1114 
1115  AliLandauGaus::EnableSigmaShift(fDoShift); // Disable sigma-shift
1116  AliLandauGausFitter f(xMin, xMax, 10);
1117  // f.SetDebug(true);
1118 
1119  TF1* r = f.FitNParticle(dist, n);
1120  if (!r) return;
1121  r->SetName("fit");
1122 
1123  // Get the status for the covariance calculation:
1124  // 0: Not calculatuted
1125  // 1: Approximate
1126  // 2: Forced positive definite
1127  // 3: Accurate
1128  Int_t covStatus =
1129  static_cast<TFitResult*>(f.GetFitResults().At(n-1))->CovMatrixStatus();
1130  if (covStatus < 2) ::Warning("", "Fit may not be good");
1131 
1132  r->SetRange(xMin, xMax);
1133  res = static_cast<TF1*>(r->Clone());
1134  res->SetLineColor(kPink+10);
1135  res->SetLineStyle(1);
1136  res->SetLineWidth(3);
1137  dist->GetListOfFunctions()->Add(res);
1138 
1139  Double_t min = 0.1*dist->GetMinimum();
1140  TF1* old = src; // Tuck away src
1141  src = res; // Use fit now
1142  for (Int_t i = 1; i <= n; i++) {
1143  Double_t rDelta = Delta();
1144  Double_t rxi = Xi();
1145  Double_t rsigma = Sigma();
1146  TF1* comp = Fi(i);
1147  comp->SetLineWidth(2);
1148  comp->SetLineStyle(1);
1149 
1150  // Modifies arguments!
1151  AliLandauGaus::IPars(i, rDelta, rxi, rsigma);
1152  TLine* l = VerticalLine(1, comp->GetLineColor(), rDelta, min,
1153  1.1*old->Eval(rDelta));
1154  dist->GetListOfFunctions()->Add(comp);
1155  dist->GetListOfFunctions()->Add(l);
1156  }
1157  src = old; // restore
1158 
1159  }
1165  void DrawInPad(TVirtualPad* p)
1166  {
1167  // Double_t scale = src->GetMaximum();
1168  // src->SetParameter(0, 1/scale);
1169  p->Clear();
1170  p->SetTopMargin(0.0);
1171  p->SetRightMargin(0.0);
1172  p->SetBottomMargin(0.0);
1173  p->SetLeftMargin(0.0);
1174  p->SetFillColor(kGray);
1175  p->Divide(2,1,0.001,0.001);
1176  TVirtualPad* q = p->cd(1);
1177  q->SetLogy();
1178  q->SetTopMargin(0.01);
1179  q->SetRightMargin(0.02);
1180  q->SetFillColor(kWhite);
1181  gStyle->SetOptStat(0);
1182  gStyle->SetOptFit(11111);
1183  gStyle->SetOptTitle(false);
1184  gStyle->SetStatY(1-q->GetTopMargin());
1185 
1186  dist->DrawCopy("hist");
1187  dist->DrawCopy("e same");
1188  AliLandauGaus::EnableSigmaShift(1); // Enable sigma-shift
1189  src->DrawClone("same");
1190  AliLandauGaus::EnableSigmaShift(fDoShift); // Disable sigma-shift
1191  if (!res) return;
1192 
1193 
1194  q = p->cd(2);
1195  q->SetTopMargin(0.01);
1196  q->SetRightMargin(0.02);
1197  q->SetFillColor(kWhite);
1198  Int_t nPar = src->GetNpar();
1199  pars = new TH1D("pars", "#chi^{2}/#nu & #Deltap_{i}/#deltap_{i}",
1200  nPar+1, 0, nPar+1);
1201  pars->SetDirectory(0);
1202  pars->SetFillColor(kCyan+2);
1203  pars->SetFillStyle(3001);
1204  pars->SetYTitle("#Deltap_{i}/#deltap_{i}");
1205 
1206  TPaveStats* stats = new TPaveStats(0,0,0,0);
1207  stats->SetX1NDC(gStyle->GetStatX()-gStyle->GetStatW());
1208  stats->SetY1NDC(gStyle->GetStatY()-gStyle->GetStatH());
1209  stats->SetX2NDC(gStyle->GetStatX());
1210  stats->SetY2NDC(gStyle->GetStatY());
1211  stats->SetTextFont(42);
1212 
1213  stats->SetBorderSize(1);
1214  stats->SetFillColor(kWhite);
1215  stats->AddText("");
1216  pars->GetListOfFunctions()->Add(stats);
1217  TString fmt(Form("%%s = %%%s", stats->GetFitFormat()));
1218  for (Int_t i = 0; i < nPar; i++) {
1219  Double_t low, high;
1220  res->GetParLimits(i, low, high);
1221  Double_t pSrc = src->GetParameter(i);
1222  Double_t pRes = res->GetParameter(i);
1223  Double_t eRes = res->GetParError(i);
1224  Bool_t fix = (low == high && ((pRes == 0 && low == 1) || low==pRes));
1225  Double_t diff = fix ? 0 : TMath::Abs(pSrc-pRes)/eRes;
1226  pars->GetXaxis()->SetBinLabel(i+1,res->GetParName(i));
1227  pars->SetBinContent(i+1, diff);
1228  stats->AddText(Form(fmt.Data(), res->GetParName(i), pSrc));
1229  }
1230  pars->GetXaxis()->SetBinLabel(nPar+1,"#chi^{2}/#nu");
1231  pars->SetBinContent(nPar+1,res->GetChisquare() / res->GetNDF());
1232  pars->DrawCopy();
1233  q->Modified();
1234  q->Update();
1235 
1236  p->Modified();
1237  p->Update();
1238  p->cd();
1239 
1240  PrintFs(false);
1241  }
1247  void PrintFs(Bool_t full=false) {
1248  if (full) {
1249  Int_t nPar = src->GetNpar();
1250  Printf("%-2s %-10s | %-8s | %8s %-9s | %9.5s",
1251  "#", "Name", "Source", "Fit", "Error", "delta");
1252  for (Int_t i = 0; i < nPar; i++) {
1253  Double_t low, high;
1254  res->GetParLimits(i, low, high);
1255  Double_t pSrc = src->GetParameter(i);
1256  Double_t pRes = res->GetParameter(i);
1257  Double_t eRes = res->GetParError(i);
1258  Bool_t fix = (low == high && ((pRes == 0 && low==1) || low==pRes));
1259  Double_t diff = fix ? 0 : TMath::Abs(pSrc-pRes)/eRes;
1260  Printf("%2d %10s | %8.4f | %8.4f %9.5f | %9.5f %s",
1261  i, src->GetParName(i), pSrc, pRes, eRes, diff,
1262  (fix ? "fixed" : (diff < 2 ? "ok" : "bad")));
1263  }
1264  Printf("chi^2/nu = %8.4f/%-3d = %9.5f",
1265  res->GetChisquare(), res->GetNDF(),
1266  res->GetChisquare() / res->GetNDF());
1267  }
1268  else
1269  Printf("%8.4f", res->GetChisquare() / res->GetNDF());
1270  }
1276  virtual void WriteOut(TDirectory* d)
1277  {
1278  d->cd();
1279  src->Write();
1280  dist->Write();
1281  if (res) res->Write();
1282  if (pars) pars->Write();
1283  }
1284 
1285  };
1286  //__________________________________________________________________
1287  TestF::Test* MakeTest() { return new Test(); }
1288  const char* Prefix() const { return "fit"; }
1289  void PreLoop(UShort_t,Int_t) {
1290  Printf("%-8s | %-8s | %-8s | %-8s ", "Delta_p", "xi", "sigma", "chi^2/nu");
1291  }
1292  void Step(UShort_t,Int_t,Int_t,Int_t,TestF::Test*,TNtuple*) {}
1293  void PostLoop(UShort_t,TCanvas* c,TNtuple*,TFile*) { PrintCanvas(c,0,0); }
1294 };
1295 #endif
1296 //
1297 // EOF
1298 //
1299 
1300 
void ScanOne(Bool_t scanSigma, Int_t n, const Double_t *values, Int_t maxN=10)
Definition: TestF.C:407
Double_t N() const
Definition: TestF.C:114
double Double_t
Definition: External.C:58
const char * title
Definition: MakeQAPdf.C:26
TCanvas * canvas
Definition: DrawAnaELoss.C:28
static TF1 * MakeFn(Double_t c, Double_t delta, Double_t xi, Double_t sigma, Double_t sigma_n, Int_t n, const Double_t *a, Double_t xmin, Double_t xmax)
const char * fmt
static TLine * VerticalLine(Int_t i, Color_t col, Double_t x, Double_t low, Double_t high)
Definition: TestF.C:195
virtual void PreLoop(UShort_t mode, Int_t nVal)=0
Double_t SigmaN() const
Definition: TestF.C:110
TSystem * gSystem
Test()
Definition: TestF.C:37
virtual ~Test()
Definition: TestF.C:46
TF1 * Fi(Double_t c, Double_t delta, Double_t xi, Double_t sigma, Double_t sigmaN, Int_t i) const
Definition: TestF.C:140
virtual void WriteOut(TDirectory *d)=0
Test(const Test &o)
Definition: TestF.C:41
TCanvas * c
Definition: TestFitELoss.C:172
virtual void DoRun()=0
void ScanOne(Bool_t scanSigma, const TArrayD &values, Int_t maxN)
Definition: TestF.C:422
virtual ~TestF()
Definition: TestF.C:203
Test * DoRealistic(TVirtualPad *pad=0)
Definition: TestF.C:286
void Scan(UShort_t mode, Int_t n, const Double_t *values, Int_t maxN)
Definition: TestF.C:459
Declaration and implementation of Landau-Gauss distributions.
TRandom * gRandom
virtual TNtuple * MakeNTuple()
Definition: TestF.C:309
Double_t Delta() const
Definition: TestF.C:98
virtual void DrawInPad(TVirtualPad *p)=0
Definition: TestF.C:25
Double_t * sigma
int Int_t
Definition: External.C:63
Double_t XMax() const
Definition: TestF.C:127
virtual const char * Prefix() const =0
Test * DoOne(const Double_t c, const Double_t delta, const Double_t xi, const Double_t sigma, const Int_t n, const Double_t *a, const Double_t xMin, const Double_t xMax, TVirtualPad *pad=0)
Definition: TestF.C:225
static TF1 * MakeFi(Double_t c, Double_t delta, Double_t xi, Double_t sigma, Double_t sigma_n, Int_t i, Double_t xmin, Double_t xmax)
Definition: External.C:212
void StepIt(UShort_t mode, Double_t xi, Double_t sigma, Int_t maxN, Int_t i, Int_t j, Int_t n, TCanvas *can, TNtuple *nt, TDirectory *p, const char *pre, Double_t v)
Definition: TestF.C:555
Double_t C() const
Definition: TestF.C:94
Int_t mode
Definition: anaM.C:40
TF1 * src
Definition: TestF.C:33
Double_t Sigma() const
Definition: TestF.C:106
TF1 * Fi(Int_t i)
Definition: TestF.C:152
static Color_t GetIColor(Int_t i)
virtual void PostLoop(UShort_t mode, TCanvas *can, TNtuple *nt, TFile *out)=0
static void IPars(Int_t i, Double_t &delta, Double_t &xi, Double_t &sigma)
Int_t nbinsx
void LoopOne(UShort_t mode, Int_t n, const Double_t *values, Int_t maxN, TCanvas *can, TNtuple *nt, TFile *out)
Definition: TestF.C:494
Double_t XMin() const
Definition: TestF.C:123
void LoopTwo(Int_t n, const Double_t *values, Int_t maxN, TCanvas *can, TNtuple *nt, TFile *out)
Definition: TestF.C:522
void ScanTwo(Int_t n, const Double_t *values, Int_t maxN=10)
Definition: TestF.C:435
void PrintCanvas(TCanvas *c, Double_t xi, Double_t sigma)
Definition: TestF.C:344
virtual void Step(UShort_t mode, Int_t i, Int_t j, Int_t n, Test *test, TNtuple *nt)=0
static Bool_t EnableSigmaShift(Short_t val=-1)
TFile * file
Test * DoUnit(const Double_t sigma, const Double_t xi=1, const Int_t n=10, TVirtualPad *pad=0)
Definition: TestF.C:260
unsigned short UShort_t
Definition: External.C:28
Double_t A(Int_t i) const
Definition: TestF.C:118
Double_t Xi() const
Definition: TestF.C:102
virtual void SetupScan(UShort_t mode, TFile *&file, TCanvas *&canvas, TNtuple *&nt)
Definition: TestF.C:318
Test & operator=(const Test &o)
Definition: TestF.C:50
bool Bool_t
Definition: External.C:53
virtual Test * MakeTest()=0
void Run(const Double_t c, const Double_t delta, const Double_t xi, const Double_t sigma, const Int_t n, const Double_t *a, const Double_t xMin, const Double_t xMax)
Definition: TestF.C:67
Definition: External.C:196
void ScanTwo(const TArrayD &values, Int_t maxN)
Definition: TestF.C:447
Declaration and implementation of fitter of Landau-Gauss distributions to energy loss spectra...