PWGLF/FORWARD/analysis2/corrs/TestF.C

   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