[root] / trunk / hist / hist / src / TMultiGraph.cxx Repository:
ViewVC logotype

Diff of /trunk/hist/hist/src/TMultiGraph.cxx

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 24702, Tue Jul 8 12:01:46 2008 UTC revision 25487, Mon Sep 22 12:44:13 2008 UTC
# Line 16  Line 16 
16  #include "TVirtualPad.h"  #include "TVirtualPad.h"
17  #include "Riostream.h"  #include "Riostream.h"
18  #include "TVirtualFitter.h"  #include "TVirtualFitter.h"
19    #include "TPluginManager.h"
20  #include "TClass.h"  #include "TClass.h"
21  #include "TMath.h"  #include "TMath.h"
22    #include "TSystem.h"
23  #include <stdlib.h>  #include <stdlib.h>
24    
25  #include <ctype.h>  #include <ctype.h>
# Line 251  Line 253 
253    
254    
255  //______________________________________________________________________________  //______________________________________________________________________________
256  Int_t TMultiGraph::Fit(TF1 *f1, Option_t *option, Option_t *, Axis_t rxmin, Axis_t rxmax)  Int_t TMultiGraph::Fit(TF1 *f1, Option_t *option, Option_t *goption, Axis_t rxmin, Axis_t rxmax)
257  {  {
258     // Fit this multigraph with function f1.     // Fit this multigraph with function f1.
259     //     //
# Line 395  Line 397 
397     //  Root > st->SetX1NDC(newx1); //new x start position     //  Root > st->SetX1NDC(newx1); //new x start position
398     //  Root > st->SetX2NDC(newx2); //new x end position     //  Root > st->SetX2NDC(newx2); //new x end position
399    
400     Int_t fitResult = 0;     return DoFit(f1,option,goption,rxmin,rxmax);  // implemented in HFitImpl.cxx
    Double_t xmin, xmax, ymin, ymax;  
    Int_t i, npar,nvpar,nparx;  
    Double_t par, we, al, bl;  
    Double_t eplus,eminus,eparab,globcc,amin,edm,errdef,werr;  
    Int_t np;  
    TF1 *fnew1;  
   
    // Check validity of function  
    if (!f1) {  
       Error("Fit", "function may not be null pointer");  
       return 0;  
    }  
    if (f1->IsZombie()) {  
       Error("Fit", "function is zombie");  
       return 0;  
    }  
   
    npar = f1->GetNpar();  
    if (npar <= 0) {  
       Error("Fit", "function %s has illegal number of parameters = %d", f1->GetName(), npar);  
       return 0;  
    }  
   
    // Check that function has same dimension as graph  
    if (f1->GetNdim() > 1) {  
       Error("Fit", "function %s is not 1-D", f1->GetName());  
       return 0;  
    }  
   
    TGraph *g;  
    TIter next(fGraphs);  
401    
    Double_t *arglist = new Double_t[100];  
    // Decode string choptin and fill fitOption structure  
    Foption_t fitOption;  
   
    TString opt = option;  
    opt.ToUpper();  
    Double_t h=0;  
    opt.ReplaceAll("ROB", "H");  
    //for robust fitting, see if # of good points is defined  
    if (opt.Contains("H=0.")) {  
       int start = opt.Index("H=0.");  
       int numpos = start + strlen("H=0.");  
       int numlen = 0;  
       int len = opt.Length();  
       while( (numpos+numlen<len) && isdigit(opt[numpos+numlen]) ) numlen++;  
       TString num = opt(numpos,numlen);  
       opt.Remove(start+strlen("H"),strlen("=0.")+numlen);  
       h = atof(num.Data());  
       h*=TMath::Power(10, -numlen);  
    }  
   
    if (opt.Contains("U")) fitOption.User    = 1;  
    if (opt.Contains("Q")) fitOption.Quiet   = 1;  
    if (opt.Contains("V")){fitOption.Verbose = 1; fitOption.Quiet   = 0;}  
    if (opt.Contains("W")) fitOption.W1      = 1;  
    if (opt.Contains("E")) fitOption.Errors  = 1;  
    if (opt.Contains("R")) fitOption.Range   = 1;  
    if (opt.Contains("N")) fitOption.Nostore = 1;  
    if (opt.Contains("0")) fitOption.Nograph = 1;  
    if (opt.Contains("+")) fitOption.Plus    = 1;  
    if (opt.Contains("B")) fitOption.Bound   = 1;  
    if (opt.Contains("C")) fitOption.Nochisq = 1;  
    if (opt.Contains("F")) fitOption.Minuit  = 1;  
    if (opt.Contains("H")) fitOption.Robust  = 1;  
   
    if (rxmax > rxmin) {  
       xmin = rxmin;  
       xmax = rxmax;  
    } else {  
       g=(TGraph *)fGraphs->First();  
       if (!g) {  
          Error("Fit", "No graphs in the multigraph");  
          return 0;  
402        }        }
       Double_t *px, *py;  
       np=g->GetN();  
       px=g->GetX();  
       py=g->GetY();  
       xmin=px[0];  
       xmax=py[np-1];  
       ymin=px[0];  
       ymax=py[np-1];  
       Double_t err0=g->GetErrorX(0);  
       Double_t errn=g->GetErrorX(np-1);  
       if (err0 > 0) xmin -= 2*err0;  
       if (errn > 0) xmax += 2*errn;  
403    
404        next.Reset();  //______________________________________________________________________________
405        while ((g = (TGraph*) next())) {  void TMultiGraph::FitPanel()
406           np=g->GetN();  {
407           px=g->GetX();  //   -*-*-*-*-*Display a panel with all histogram fit options*-*-*-*-*-*
408           py=g->GetY();  //             ==============================================
409           for (i=0; i<np; i++) {  //
410              if (px[i] < xmin) xmin = px[i];  //      See class TFitPanel for example
             if (px[i] > xmax) xmax = px[i];  
             if (py[i] < ymin) ymin = py[i];  
             if (py[i] > ymax) ymax = py[i];  
          }  
       }  
    }  
   
    ///////////////  
    //set the fitter  
    //////////////  
   
    Int_t special=f1->GetNumber();  
    Bool_t linear = f1->IsLinear();  
    if (special==299+npar)  
       linear=kTRUE;  
    if (fitOption.Bound || fitOption.User || fitOption.Errors || fitOption.Minuit)  
       linear = kFALSE;  
   
    char l[]="TLinearFitter";  
    Int_t strdiff = 0;  
    Bool_t isSet = kFALSE;  
    if (TVirtualFitter::GetFitter()){  
       //Is a fitter already set? Is it linear?  
       isSet = kTRUE;  
       strdiff = strcmp(TVirtualFitter::GetFitter()->IsA()->GetName(), l);  
    }  
    if (linear){  
       TClass *cl = TClass::GetClass("TLinearFitter");  
       if (isSet && strdiff!=0) {  
          delete TVirtualFitter::GetFitter();  
          isSet=kFALSE;  
       }  
       if (!isSet && cl) {  
          TVirtualFitter::SetFitter((TVirtualFitter *)cl->New());  
       }  
    } else {  
       if (isSet && strdiff==0){  
          delete TVirtualFitter::GetFitter();  
          isSet=kFALSE;  
       }  
       if (!isSet)  
          TVirtualFitter::SetFitter(0);  
    }  
   
    TVirtualFitter *grFitter = TVirtualFitter::Fitter(this, f1->GetNpar());  
    grFitter->Clear();  
   
    // Get pointer to the function by searching in the list of functions in ROOT  
    grFitter->SetUserFunc(f1);  
    grFitter->SetFitOption(fitOption);  
   
    // Is a Fit range specified?  
    if (fitOption.Range) {  
       f1->GetRange(xmin, xmax);  
    } else {  
       f1->SetRange(xmin, xmax);  
    }  
   
    if (linear){  
       if (fitOption.Robust)  
          fitResult = grFitter->ExecuteCommand("FitMultiGraph", &h, 0);  
       else  
          fitResult = grFitter->ExecuteCommand("FitMultiGraph", 0, 0);  
    } else {  
   
       //Int_t special = f1->GetNumber();  
       if (fitOption.Bound) special = 0;  
       if      (special == 100)      InitGaus(xmin,xmax);  
       else if (special == 400)      InitGaus(xmin,xmax);  
       else if (special == 200)      InitExpo(xmin,xmax);  
       else if (special == 299+npar) InitPolynom(xmin,xmax);  
   
       //*-*- Some initialisations  
       if (!fitOption.Verbose) {  
       arglist[0] = -1;  
       grFitter->ExecuteCommand("SET PRINT", arglist,1);  
       arglist[0] = 0;  
       grFitter->ExecuteCommand("SET NOW",   arglist,0);  
       }  
   
       /////////////////////////////////////////////////////////  
       //*-*- Set error criterion for chisquare  
       arglist[0] = TVirtualFitter::GetErrorDef();  
       if (!fitOption.User) grFitter->SetFitMethod("MultiGraphFitChisquare");  
   
   
       fitResult = grFitter->ExecuteCommand("SET ERR",arglist,1);  
       if (fitResult != 0) {  
          //   Abnormal termination, MIGRAD might not have converged on a  
          //   minimum.  
          if (!fitOption.Quiet) {  
             Warning("Fit","Abnormal termination of minimization.");  
          }  
          delete [] arglist;  
          return fitResult;  
       }  
   
       //*-*- Transfer names and initial values of parameters to Minuit  
       Int_t nfixed = 0;  
       for (i=0;i<npar;i++) {  
          par = f1->GetParameter(i);  
          f1->GetParLimits(i,al,bl);  
          if (al*bl != 0 && al >= bl) {  
             al = bl = 0;  
             arglist[nfixed] = i+1;  
             nfixed++;  
          }  
          we  = 0.3*TMath::Abs(par);  
          if (we <= TMath::Abs(par)*1e-6) we = 1;  
          grFitter->SetParameter(i,f1->GetParName(i),par,we,al,bl);  
       }  
       if(nfixed > 0)grFitter->ExecuteCommand("FIX",arglist,nfixed); // Otto  
   
       //*-*- Reset Print level  
       if (!fitOption.Quiet) {  
          if (fitOption.Verbose) { arglist[0] = 2; grFitter->ExecuteCommand("SET PRINT", arglist,1); }  
          else                   { arglist[0] = 0; grFitter->ExecuteCommand("SET PRINT", arglist,1); }  
       }  
       //*-*- Compute sum of squares of errors in the bin range  
       Bool_t hasErrors = kFALSE;  
       Double_t ex, ey, sumw2=0;  
       next.Reset();  
       while ((g = (TGraph*) next())) {  
          np=g->GetN();  
          for (i=0; i<np; i++){  
             ex=g->GetErrorX(i);  
             ey=g->GetErrorY(i);  
             if (ex > 0 || ey > 0) hasErrors=kTRUE;  
             sumw2+=ey*ey;  
          }  
       }  
   
       //*-*- Perform minimization  
   
       arglist[0] = TVirtualFitter::GetMaxIterations();  
       arglist[1] = sumw2*TVirtualFitter::GetPrecision();  
       grFitter->ExecuteCommand("MIGRAD",arglist,2);  
       if (fitOption.Errors) {  
          grFitter->ExecuteCommand("HESSE",arglist,0);  
          grFitter->ExecuteCommand("MINOS",arglist,0);  
       }  
   
       grFitter->GetStats(amin,edm,errdef,nvpar,nparx);  
       f1->SetChisquare(amin);  
       Int_t ndf = f1->GetNumberFitPoints()-npar+nfixed;  
       f1->SetNDF(ndf);  
   
       //*-*- Get return status  
       char parName[50];  
       for (i=0;i<npar;i++) {  
          grFitter->GetParameter(i,parName, par,we,al,bl);  
          if (!fitOption.Errors) werr = we;  
          else {  
             grFitter->GetErrors(i,eplus,eminus,eparab,globcc);  
             if (eplus > 0 && eminus < 0) werr = 0.5*(eplus-eminus);  
             else                         werr = we;  
          }  
          if (!hasErrors && ndf > 1) werr *= TMath::Sqrt(amin/(ndf-1));  
          f1->SetParameter(i,par);  
          f1->SetParError(i,werr);  
       }  
    }  
    // Print final values of parameters.  
    if (!fitOption.Quiet) {  
       if (fitOption.Errors) grFitter->PrintResults(4,amin);  
       else                  grFitter->PrintResults(3,amin);  
    }  
    delete [] arglist;  
411    
412     // Store fitted function in histogram functions list and draw     if (!gPad)
413          gROOT->MakeDefCanvas();
414    
415     if (!fitOption.Nostore) {     if (!gPad) {
416        if (!fFunctions) fFunctions = new TList;        Error("FitPanel", "Unable to create a default canvas");
417        if (!fitOption.Plus) {        return;
          TIter next2(fFunctions, kIterBackward);  
          TObject *obj;  
          while ((obj = next2())) {  
             if (obj->InheritsFrom(TF1::Class())){  
                obj = fFunctions->Remove(obj);  
                delete obj;  
             }  
          }  
418        }        }
       fnew1 = new TF1();  
       f1->Copy(*fnew1);  
       fFunctions->Add(fnew1);  
       fnew1->SetParent(this);  
       fnew1->Save(xmin,xmax,0,0,0,0);  
       if (fitOption.Nograph) fnew1->SetBit(TF1::kNotDraw);  
       fnew1->SetBit(TFormula::kNotGlobal);  
419    
420        if (TestBit(kCanDelete)) return fitResult;     // use plugin manager to create instance of TFitEditor
421        if (gPad) gPad->Modified();     TPluginHandler *handler = gROOT->GetPluginManager()->FindHandler("TFitEditor");
422       if (handler && handler->LoadPlugin() != -1) {
423          if (handler->ExecPlugin(2, gPad, this) == 0)
424             Error("FitPanel", "Unable to crate the FitPanel");
425     }     }
426       else
427     return fitResult;           Error("FitPanel", "Unable to find the FitPanel plug-in");
428  }  }
429    
   
430  //______________________________________________________________________________  //______________________________________________________________________________
431  Option_t *TMultiGraph::GetGraphDrawOption(const TGraph *gr) const  Option_t *TMultiGraph::GetGraphDrawOption(const TGraph *gr) const
432  {  {
# Line 970  Line 706 
706     return (TF1*)fFunctions->FindObject(name);     return (TF1*)fFunctions->FindObject(name);
707  }  }
708    
709    //______________________________________________________________________________
710    TList *TMultiGraph::GetListOfFunctions()
711    {
712       // Return pointer to list of functions
713       // if pointer is null create the list
714    
715       if (!fFunctions) fFunctions = new TList();
716       return fFunctions;
717    }
718    
719    
720  //______________________________________________________________________________  //______________________________________________________________________________
721  TAxis *TMultiGraph::GetXaxis() const  TAxis *TMultiGraph::GetXaxis() const
Legend:
Removed from v.24702  
changed lines
  Added in v.25487
Subversion Admin
 ViewVC Help
Powered by ViewVC 1.0.9