TBackCompFitter.cxx

Go to the documentation of this file.
00001 #include "TROOT.h"
00002 #include "TBackCompFitter.h"
00003 
00004 
00005 #include "TMethodCall.h"
00006 #include "TInterpreter.h"
00007 
00008 #include "Math/Util.h"
00009 
00010 #include <iostream>
00011 #include <cassert>
00012 
00013 //needed by GetCondifenceLevel
00014 #include "Math/IParamFunction.h"
00015 #include "TH1.h"
00016 #include "TH2.h"
00017 #include "TH3.h"
00018 #include "TMath.h"
00019 #include "TGraph.h"
00020 #include "TGraphErrors.h"
00021 #include "TGraph2DErrors.h"
00022 #include "TMultiGraph.h"
00023 #include "HFitInterface.h"
00024 #include "Math/Minimizer.h"
00025 #include "Fit/BinData.h"
00026 #include "Fit/UnBinData.h"
00027 #include "Fit/PoissonLikelihoodFCN.h"
00028 #include "Fit/LogLikelihoodFCN.h"
00029 #include "Fit/Chi2FCN.h"
00030 #include "Fit/FcnAdapter.h"
00031 #include "TFitResult.h"
00032 
00033 //#define DEBUG 1
00034 
00035 
00036 //______________________________________________________________________________
00037 /**
00038    Backward compatible implementation of TVirtualFitter using the the class ROOT::Fit::Fitter.
00039    This class is created after fitting an histogram (TH1), TGraph or TTree and provides in addition to the 
00040    methods of the TVirtualFitter hooks to access the fit result class (ROOT::Fit::FitResult), the fit configuration
00041    (ROOT::Fit::FitConfig) or the fit data (ROOT::Fit::FitData) using
00042    <pre>
00043    TBackCompFitter * fitter = (TBackCompFitter *) TVirtualFitter::GetFitter();
00044    ROOT::Fit::FitResult & result = fitter->GetFitResult();
00045    result.Print(std::cout);
00046    </pre>
00047    
00048    Methods for getting the confidence level or contours are also provided. 
00049    Note that after a new calls to TH1::Fit (or similar) the class will be deleted and all reference to the FitResult, FitConfig 
00050    or minimizer will be invalid. One could eventually copying  the class before issuing a new fit to avoid deleting this information
00051 */
00052 
00053 
00054 
00055 ClassImp(TBackCompFitter);
00056 
00057 
00058 
00059 TBackCompFitter::TBackCompFitter( ) : 
00060    fMinimizer(0), 
00061    fObjFunc(0),
00062    fModelFunc(0)
00063 {
00064    // Constructur needed by TVirtualFitter interface. Same behavior as default constructor.
00065    // initialize setting name and the global pointer
00066    SetName("BCFitter");
00067 }
00068 
00069 TBackCompFitter::TBackCompFitter(std::auto_ptr<ROOT::Fit::Fitter> fitter, std::auto_ptr<ROOT::Fit::FitData>  data) : 
00070    fFitData(data),
00071    fFitter(fitter),
00072    fMinimizer(0),
00073    fObjFunc(0),
00074    fModelFunc(0)
00075 {
00076    // constructor used after having fit using directly ROOT::Fit::Fitter
00077    // will create a dummy fitter copying configuration and parameter settings
00078    SetName("LastFitter");
00079 }
00080 
00081 
00082 
00083 
00084 TBackCompFitter::~TBackCompFitter() { 
00085    // data are own here
00086    //if (fFitData) delete fFitData; 
00087 
00088    if (fMinimizer) delete fMinimizer; 
00089    if (fObjFunc) delete fObjFunc; 
00090    if (fModelFunc) delete fModelFunc;
00091 }
00092 
00093 Double_t TBackCompFitter::Chisquare(Int_t npar, Double_t *params) const {
00094    // do chisquare calculations in case of likelihood fits 
00095    // do evaluation a the minimum only 
00096    const std::vector<double> & minpar = fFitter->Result().Parameters(); 
00097    assert (npar == (int) minpar.size() ); 
00098    double diff = 0; 
00099    double s = 0; 
00100    for (int i =0; i < npar; ++i) { 
00101       diff += std::abs( params[i] - minpar[i] );  
00102       s += minpar[i]; 
00103    }
00104 
00105    if (diff > s * 1.E-12 ) Warning("Chisquare","given parameter values are not at minimum - chi2 at minimum is returned"); 
00106    return fFitter->Result().Chi2(); 
00107 }
00108 
00109 void TBackCompFitter::Clear(Option_t*) {
00110    // clear resources for consecutive fits
00111 
00112    // need to do something here ??? to be seen
00113    
00114    
00115 }
00116 
00117 
00118 
00119 
00120 
00121 Int_t TBackCompFitter::ExecuteCommand(const char *command, Double_t *args, Int_t nargs){        
00122    // execute the command (Fortran Minuit compatible interface)
00123    
00124 #ifdef DEBUG
00125    std::cout<<"Execute command= "<<command<<std::endl;
00126 #endif
00127 
00128    int nfcn = GetMaxIterations();  
00129    double edmval = GetPrecision(); 
00130 
00131    // set also number of parameters in obj function
00132    DoSetDimension(); 
00133 
00134    TString scommand(command); 
00135    scommand.ToUpper();
00136 
00137    // MIGRAD 
00138    if (scommand.Contains("MIG")) {
00139       if (nargs > 0) nfcn = int ( args[0] );
00140       if (nargs > 1) edmval = args[1];
00141       if (!fObjFunc) { 
00142          Error("ExecuteCommand","FCN must set before executing this command"); 
00143          return -1; 
00144       }
00145 
00146       fFitter->Config().SetMinimizer(GetDefaultFitter(), "Migrad");
00147       bool ret = fFitter->FitFCN(*fObjFunc); 
00148       return  (ret) ? 0 : -1;
00149       
00150       
00151       
00152    } 
00153    //Minimize
00154    if (scommand.Contains("MINI")) {
00155 
00156       if (nargs > 0) nfcn = int ( args[0] );
00157       if (nargs > 1) edmval = args[1];
00158 
00159       fFitter->Config().SetMinimizer(GetDefaultFitter(), "Minimize");
00160       if (!fObjFunc) { 
00161          Error("ExecuteCommand","FCN must set before executing this command"); 
00162          return -1; 
00163       }
00164       bool ret = fFitter->FitFCN(*fObjFunc); 
00165       return  (ret) ? 0 : -1;
00166    }
00167    //Simplex
00168    if (scommand.Contains("SIM")) {
00169       
00170       if (nargs > 0) nfcn = int ( args[0] );
00171       if (nargs > 1) edmval = args[1];
00172       if (!fObjFunc) { 
00173          Error("ExecuteCommand","FCN must set before executing this command"); 
00174          return -1; 
00175       }
00176 
00177       fFitter->Config().SetMinimizer(GetDefaultFitter(), "Simplex");
00178       bool ret = fFitter->FitFCN(*fObjFunc); 
00179       return  (ret) ? 0 : -1;
00180    }
00181    //SCan
00182    if (scommand.Contains("SCA")) {
00183       
00184       if (nargs > 0) nfcn = int ( args[0] );
00185       if (nargs > 1) edmval = args[1];
00186       if (!fObjFunc) { 
00187          Error("ExecuteCommand","FCN must set before executing this command"); 
00188          return -1; 
00189       }
00190 
00191       fFitter->Config().SetMinimizer(GetDefaultFitter(), "Scan");
00192       bool ret = fFitter->FitFCN(*fObjFunc); 
00193       return  (ret) ? 0 : -1;
00194    }
00195    // MINOS 
00196    else if (scommand.Contains("MINO"))   {
00197 
00198       if (fFitter->Config().MinosErrors() ) return 0; 
00199 
00200       if (!fObjFunc) { 
00201          Error("ExecuteCommand","FCN must set before executing this command"); 
00202          return -1; 
00203       }
00204       // do only MINOS. need access to minimizer. For the moment re-run fitting with minos options 
00205       fFitter->Config().SetMinosErrors(true);
00206       // set new parameter values
00207 
00208       fFitter->Config().SetMinimizer(GetDefaultFitter(), "Migrad"); // redo -minimization with Minos
00209       bool ret = fFitter->FitFCN(*fObjFunc); 
00210       return  (ret) ? 0 : -1;
00211 
00212    } 
00213    //HESSE
00214    else if (scommand.Contains("HES"))   {
00215 
00216       if (fFitter->Config().ParabErrors() ) return 0; 
00217 
00218       if (!fObjFunc) { 
00219          Error("ExecuteCommand","FCN must set before executing this command"); 
00220          return -1; 
00221       }
00222 
00223       // do only HESSE. need access to minimizer. For the moment re-run fitting with hesse options 
00224       fFitter->Config().SetParabErrors(true);
00225       fFitter->Config().SetMinimizer(GetDefaultFitter(), "Migrad"); // redo -minimization with Minos
00226       bool ret = fFitter->FitFCN(*fObjFunc); 
00227       return  (ret) ? 0 : -1;
00228    } 
00229    
00230    // FIX 
00231    else if (scommand.Contains("FIX"))   {
00232       for(int i = 0; i < nargs; i++) {
00233          FixParameter(int(args[i])-1);
00234       }
00235       return 0;
00236    } 
00237    // SET LIMIT (upper and lower)
00238    else if (scommand.Contains("SET LIM"))   {
00239       if (nargs < 3) { 
00240          Error("ExecuteCommand","Invalid parameters given in SET LIMIT");
00241          return -1; 
00242       }
00243       int ipar = int(args[0]);
00244       if (!ValidParameterIndex(ipar) )  return -1;   
00245       double low = args[1];
00246       double up = args[2];
00247       fFitter->Config().ParSettings(ipar).SetLimits(low,up);
00248       return 0; 
00249    } 
00250    // SET PRINT
00251    else if (scommand.Contains("SET PRIN"))   {
00252       if (nargs < 1) return -1;  
00253       fFitter->Config().MinimizerOptions().SetPrintLevel(int(args[0]) );
00254       return 0; 
00255    } 
00256    // SET ERR
00257    else if (scommand.Contains("SET ERR"))   {
00258       if (nargs < 1) return -1;  
00259       fFitter->Config().MinimizerOptions().SetPrintLevel(int( args[0]) );
00260       return 0; 
00261    } 
00262    // SET STRATEGY
00263    else if (scommand.Contains("SET STR"))   {
00264       if (nargs < 1) return -1;  
00265       fFitter->Config().MinimizerOptions().SetStrategy(int(args[0]) );
00266       return 0; 
00267    } 
00268    //SET GRAD (not impl.) 
00269    else if (scommand.Contains("SET GRA"))   {
00270       //     not yet available
00271       //     fGradient = true;
00272       return -1;
00273    } 
00274    //SET NOW (not impl.) 
00275    else if (scommand.Contains("SET NOW"))   {
00276       //     no warning (works only for TMinuit)
00277       //     fGradient = true;
00278       return -1;
00279    } 
00280    // CALL FCN
00281    else if (scommand.Contains("CALL FCN"))   {
00282       //     call fcn function (global pointer to free function)
00283 
00284       if (nargs < 1 || fFCN == 0 ) return -1;
00285       int npar = fObjFunc->NDim();
00286       // use values in fit result if existing  otherwise in ParameterSettings
00287       std::vector<double> params(npar); 
00288       for (int i = 0; i < npar; ++i) 
00289          params[i] = GetParameter(i); 
00290 
00291       double fval = 0;
00292       (*fFCN)(npar, 0, fval, &params[0],int(args[0]) ) ;
00293       return 0; 
00294    } 
00295    else {
00296       // other commands passed 
00297       Error("ExecuteCommand","Invalid or not supported command given %s",command);
00298       return -1;
00299    }
00300    
00301    
00302 }
00303 
00304 bool TBackCompFitter::ValidParameterIndex(int ipar) const  { 
00305    // check if ipar is a valid parameter index
00306    int nps  = fFitter->Config().ParamsSettings().size(); 
00307    if (ipar  < 0 || ipar >= nps ) { 
00308       std::string msg = ROOT::Math::Util::ToString(ipar) + " is an invalid Parameter index";
00309       Error("ValidParameterIndex","%s",msg.c_str());
00310       return false;
00311    } 
00312    return true; 
00313 }
00314 
00315 void TBackCompFitter::FixParameter(Int_t ipar) {
00316    // fix the paramter
00317    //   std::cout<<"FixParameter"<<std::endl;
00318    if (ValidParameterIndex(ipar) )    
00319       fFitter->Config().ParSettings(ipar).Fix();
00320 }
00321 
00322 
00323 
00324 void TBackCompFitter::GetConfidenceIntervals(Int_t n, Int_t ndim, const Double_t *x, Double_t *ci, Double_t cl)
00325 {
00326 //Computes point-by-point confidence intervals for the fitted function
00327 //Parameters:
00328 //n - number of points
00329 //ndim - dimensions of points
00330 //x - points, at which to compute the intervals, for ndim > 1 
00331 //    should be in order: (x0,y0, x1, y1, ... xn, yn)
00332 //ci - computed intervals are returned in this array
00333 //cl - confidence level, default=0.95
00334 //NOTE, that the intervals are approximate for nonlinear(in parameters) models
00335    
00336    if (!fFitter->Result().IsValid()) { 
00337       Error("GetConfidenceIntervals","Cannot compute confidence intervals with an invalide fit result");
00338       return; 
00339    }
00340    
00341    fFitter->Result().GetConfidenceIntervals(n,ndim,1,x,ci,cl);         
00342 }
00343 
00344 void TBackCompFitter::GetConfidenceIntervals(TObject *obj, Double_t cl)
00345 {
00346 //Computes confidence intervals at level cl. Default is 0.95
00347 //The TObject parameter can be a TGraphErrors, a TGraph2DErrors or a TH1,2,3.
00348 //For Graphs, confidence intervals are computed for each point,
00349 //the value of the graph at that point is set to the function value at that
00350 //point, and the graph y-errors (or z-errors) are set to the value of
00351 //the confidence interval at that point.
00352 //For Histograms, confidence intervals are computed for each bin center
00353 //The bin content of this bin is then set to the function value at the bin
00354 //center, and the bin error is set to the confidence interval value.
00355 //NOTE: confidence intervals are approximate for nonlinear models!
00356 //
00357 //Allowed combinations:
00358 //Fitted object               Passed object
00359 //TGraph                      TGraphErrors, TH1
00360 //TGraphErrors, AsymmErrors   TGraphErrors, TH1
00361 //TH1                         TGraphErrors, TH1
00362 //TGraph2D                    TGraph2DErrors, TH2
00363 //TGraph2DErrors              TGraph2DErrors, TH2
00364 //TH2                         TGraph2DErrors, TH2
00365 //TH3                         TH3
00366 
00367    if (!fFitter->Result().IsValid() ) { 
00368       Error("GetConfidenceIntervals","Cannot compute confidence intervals with an invalide fit result");
00369       return; 
00370    }
00371 
00372    // get data dimension from fit object
00373    int datadim = 1; 
00374    TObject * fitobj = GetObjectFit(); 
00375    if (!fitobj) { 
00376       Error("GetConfidenceIntervals","Cannot compute confidence intervals without a fitting object");
00377       return; 
00378    }
00379 
00380    if (fitobj->InheritsFrom(TGraph2D::Class())) datadim = 2; 
00381    if (fitobj->InheritsFrom(TH1::Class())) { 
00382       TH1 * h1 = dynamic_cast<TH1*>(fitobj); 
00383       assert(h1 != 0); 
00384       datadim = h1->GetDimension(); 
00385    } 
00386 
00387    if (datadim == 1) { 
00388       if (!obj->InheritsFrom(TGraphErrors::Class()) && !obj->InheritsFrom(TH1::Class() ) )  {
00389          Error("GetConfidenceIntervals", "Invalid object passed for storing confidence level data, must be a TGraphErrors or a TH1");
00390          return; 
00391       }
00392    } 
00393    if (datadim == 2) { 
00394       if (!obj->InheritsFrom(TGraph2DErrors::Class()) && !obj->InheritsFrom(TH2::Class() ) )  {
00395          Error("GetConfidenceIntervals", "Invalid object passed for storing confidence level data, must be a TGraph2DErrors or a TH2");
00396          return; 
00397       }
00398    }
00399    if (datadim == 3) { 
00400       if (!obj->InheritsFrom(TH3::Class() ) )  {
00401          Error("GetConfidenceIntervals", "Invalid object passed for storing confidence level data, must be a TH3");
00402          return; 
00403       }
00404    }
00405 
00406    // fill bin data (for the moment use all ranges) according to object passed
00407    ROOT::Fit::BinData data; 
00408    data.Opt().fUseEmpty = true; // need to use all bins of given histograms
00409    // call appropriate function according to type of object
00410    if (obj->InheritsFrom(TGraph::Class()) ) 
00411       ROOT::Fit::FillData(data, dynamic_cast<TGraph *>(obj) ); 
00412    else if (obj->InheritsFrom(TGraph2D::Class()) ) 
00413       ROOT::Fit::FillData(data, dynamic_cast<TGraph2D *>(obj) ); 
00414 //    else if (obj->InheritsFrom(TMultiGraph::Class()) ) 
00415 //       ROOT::Fit::FillData(data, dynamic_cast<TMultiGraph *>(obj) ); 
00416    else if (obj->InheritsFrom(TH1::Class()) ) 
00417       ROOT::Fit::FillData(data, dynamic_cast<TH1 *>(obj) ); 
00418    
00419 
00420    unsigned int n = data.Size(); 
00421 
00422    std::vector<double> ci( n ); 
00423 
00424    fFitter->Result().GetConfidenceIntervals(data,&ci[0],cl);         
00425 
00426    const ROOT::Math::IParamMultiFunction * func =  fFitter->Result().FittedFunction(); 
00427    assert(func != 0); 
00428 
00429    // fill now the object with cl data
00430    for (unsigned int i = 0; i < n; ++i) {
00431       const double * x = data.Coords(i); 
00432       double y = (*func)( x ); // function is evaluated using its  parameters
00433 
00434       if (obj->InheritsFrom(TGraphErrors::Class()) ) { 
00435          TGraphErrors * gr = dynamic_cast<TGraphErrors *> (obj); 
00436          assert(gr != 0); 
00437          gr->SetPoint(i, *x, y); 
00438          gr->SetPointError(i, 0, ci[i]); 
00439       }
00440       if (obj->InheritsFrom(TGraph2DErrors::Class()) ) { 
00441          TGraph2DErrors * gr = dynamic_cast<TGraph2DErrors *> (obj); 
00442          assert(gr != 0); 
00443          gr->SetPoint(i, x[0], x[1], y); 
00444          gr->SetPointError(i, 0, 0, ci[i]); 
00445       }
00446       if (obj->InheritsFrom(TH1::Class()) ) { 
00447          TH1 * h1 = dynamic_cast<TH1 *> (obj); 
00448          assert(h1 != 0); 
00449          int ibin = 0; 
00450          if (datadim == 1) ibin = h1->FindBin(*x); 
00451          if (datadim == 2) ibin = h1->FindBin(x[0],x[1]); 
00452          if (datadim == 3) ibin = h1->FindBin(x[0],x[1],x[2]); 
00453          h1->SetBinContent(ibin, y); 
00454          h1->SetBinError(ibin, ci[i]); 
00455       }
00456    }
00457 
00458 }
00459 
00460 Double_t* TBackCompFitter::GetCovarianceMatrix() const {
00461    // get the error matrix in a pointer to a NxN array.  
00462    // excluding the fixed parameters 
00463 
00464    unsigned int nfreepar =   GetNumberFreeParameters();
00465    unsigned int ntotpar =   GetNumberTotalParameters();
00466    
00467    if (fCovar.size() !=  nfreepar*nfreepar ) 
00468       fCovar.resize(nfreepar*nfreepar);
00469 
00470    if (!fFitter->Result().IsValid() ) { 
00471       Warning("GetCovarianceMatrix","Invalid fit result");
00472       return 0; 
00473    }
00474 
00475    unsigned int l = 0; 
00476    for (unsigned int i = 0; i < ntotpar; ++i) { 
00477       if (fFitter->Config().ParSettings(i).IsFixed() ) continue;
00478       unsigned int m = 0; 
00479       for (unsigned int j = 0; j < ntotpar; ++j) {
00480          if (fFitter->Config().ParSettings(j).IsFixed() ) continue;
00481          unsigned int index = nfreepar*l + m;
00482          assert(index < fCovar.size() );
00483          fCovar[index] = fFitter->Result().CovMatrix(i,j);
00484          m++;
00485       }
00486       l++;
00487    }
00488    return &(fCovar.front());
00489 }
00490 
00491 Double_t TBackCompFitter::GetCovarianceMatrixElement(Int_t i, Int_t j) const {
00492    // get error matrix element (return all zero if matrix is not available)
00493 
00494    unsigned int np2 = fCovar.size();
00495    unsigned int npar = GetNumberFreeParameters(); 
00496    if ( np2 == 0 || np2 != npar *npar ) { 
00497       double * c = GetCovarianceMatrix();
00498       if (c == 0) return 0;  
00499    }
00500    return fCovar[i*npar + j];  
00501 }
00502 
00503 
00504 Int_t TBackCompFitter::GetErrors(Int_t ipar,Double_t &eplus, Double_t &eminus, Double_t &eparab, Double_t &globcc) const {
00505    // get fit errors 
00506 
00507    if (!ValidParameterIndex(ipar) )   return -1; 
00508    
00509    const ROOT::Fit::FitResult & result = fFitter->Result(); 
00510    if (!result.IsValid() ) { 
00511       Warning("GetErrors","Invalid fit result");
00512       return -1; 
00513    }
00514 
00515    eparab = result.Error(ipar); 
00516    eplus = result.UpperError(ipar); 
00517    eminus = result.LowerError(ipar); 
00518    globcc = result.GlobalCC(ipar); 
00519    return 0;
00520 }
00521 
00522 Int_t TBackCompFitter::GetNumberTotalParameters() const {
00523    // number of total parameters 
00524    return fFitter->Result().NTotalParameters();  
00525 }
00526 Int_t TBackCompFitter::GetNumberFreeParameters() const {
00527    // number of variable parameters
00528    return fFitter->Result().NFreeParameters();  
00529 }
00530 
00531 
00532 Double_t TBackCompFitter::GetParError(Int_t ipar) const {
00533    // parameter error
00534    if (fFitter->Result().IsEmpty() ) {
00535       if (ValidParameterIndex(ipar) )  return  fFitter->Config().ParSettings(ipar).StepSize();
00536       else return 0; 
00537    }
00538    return fFitter->Result().Error(ipar);  
00539 }
00540 
00541 Double_t TBackCompFitter::GetParameter(Int_t ipar) const {
00542    // parameter value
00543    if (fFitter->Result().IsEmpty() ) {
00544       if (ValidParameterIndex(ipar) )  return  fFitter->Config().ParSettings(ipar).Value();
00545       else return 0; 
00546    }
00547    return fFitter->Result().Value(ipar);  
00548 }
00549 
00550 Int_t TBackCompFitter::GetParameter(Int_t ipar,char *name,Double_t &value,Double_t &verr,Double_t &vlow, Double_t &vhigh) const {
00551    // get all parameter info (name, value, errors) 
00552    if (!ValidParameterIndex(ipar) )    {
00553       return -1; 
00554    }
00555    const std::string & pname = fFitter->Config().ParSettings(ipar).Name(); 
00556    const char * c = pname.c_str(); 
00557    std::copy(c,c + pname.size(),name);
00558 
00559    if (fFitter->Result().IsEmpty() ) { 
00560       value = fFitter->Config().ParSettings(ipar).Value(); 
00561       verr  = fFitter->Config().ParSettings(ipar).Value();  // error is step size in this case 
00562       vlow  = fFitter->Config().ParSettings(ipar).LowerLimit();  // vlow is lower limit in this case 
00563       vhigh   = fFitter->Config().ParSettings(ipar).UpperLimit();  // vlow is lower limit in this case 
00564       return 1; 
00565    }
00566    else { 
00567       value =  fFitter->Result().Value(ipar);  
00568       verr = fFitter->Result().Error(ipar);  
00569       vlow = fFitter->Result().LowerError(ipar);  
00570       vhigh = fFitter->Result().UpperError(ipar);  
00571    }
00572    return 0; 
00573 }
00574 
00575 const char *TBackCompFitter::GetParName(Int_t ipar) const {
00576    //   return name of parameter ipar
00577    if (!ValidParameterIndex(ipar) )    {
00578       return 0; 
00579    }
00580    return fFitter->Config().ParSettings(ipar).Name().c_str(); 
00581 }
00582 
00583 Int_t TBackCompFitter::GetStats(Double_t &amin, Double_t &edm, Double_t &errdef, Int_t &nvpar, Int_t &nparx) const {
00584    // get fit statistical information
00585    const ROOT::Fit::FitResult & result = fFitter->Result(); 
00586    amin = result.MinFcnValue(); 
00587    edm = result.Edm(); 
00588    errdef = fFitter->Config().MinimizerOptions().ErrorDef(); 
00589    nvpar = result.NFreeParameters();  
00590    nparx = result.NTotalParameters();  
00591    return 0;
00592 }
00593 
00594 Double_t TBackCompFitter::GetSumLog(Int_t) {
00595    //   sum of log . Un-needed
00596    Warning("GetSumLog","Dummy  method - returned 0"); 
00597    return 0.;
00598 }
00599 
00600 
00601 Bool_t TBackCompFitter::IsFixed(Int_t ipar) const {
00602    // query if parameter ipar is fixed
00603    if (!ValidParameterIndex(ipar) )    {
00604       return false; 
00605    }
00606    return fFitter->Config().ParSettings(ipar).IsFixed(); 
00607 }
00608 
00609 
00610 void TBackCompFitter::PrintResults(Int_t level, Double_t ) const {
00611    // print the fit result
00612    // use PrintResults function in case of Minuit for old -style printing
00613    if (fFitter->GetMinimizer() && fFitter->Config().MinimizerType() == "Minuit")
00614       fFitter->GetMinimizer()->PrintResults();
00615    else { 
00616       if (level > 0) fFitter->Result().Print(std::cout); 
00617       if (level > 3)  fFitter->Result().PrintCovMatrix(std::cout);    
00618    }
00619    // need to print minos errors and globalCC + other info
00620 }
00621 
00622 void TBackCompFitter::ReleaseParameter(Int_t ipar) {
00623    // release a fit parameter
00624    if (ValidParameterIndex(ipar) )    
00625       fFitter->Config().ParSettings(ipar).Release(); 
00626 }
00627 
00628 
00629 
00630 void TBackCompFitter::SetFitMethod(const char *) {
00631    // set fit method (i.e. chi2 or likelihood)
00632    // according to the method the appropriate FCN function will be created   
00633    Info("SetFitMethod","non supported method");
00634 }
00635 
00636 Int_t TBackCompFitter::SetParameter(Int_t ipar,const char *parname,Double_t value,Double_t verr,Double_t vlow, Double_t vhigh) {
00637    // set (add) a new fit parameter passing initial value,  step size (verr) and parametr limits
00638    // if vlow > vhigh the parameter is unbounded
00639    // if the stepsize (verr) == 0 the parameter is treated as fixed   
00640 
00641    std::vector<ROOT::Fit::ParameterSettings> & parlist = fFitter->Config().ParamsSettings(); 
00642    if ( ipar >= (int) parlist.size() ) parlist.resize(ipar+1); 
00643    ROOT::Fit::ParameterSettings ps(parname, value, verr); 
00644    if (verr == 0) ps.Fix(); 
00645    if (vlow < vhigh) ps.SetLimits(vlow, vhigh); 
00646    parlist[ipar] = ps; 
00647    return 0; 
00648 }
00649 
00650 // static method evaluating FCN
00651 // void TBackCompFitter::FCN( int &, double * , double & f, double * x , int /* iflag */) { 
00652 //    // get static instance of fitter
00653 //    TBackCompFitter * fitter = dynamic_cast<TBackCompFitter *>(TVirtualFitter::GetFitter()); 
00654 //    assert(fitter); 
00655 //    if (fitter->fObjFunc == 0) fitter->RecreateFCN(); 
00656 //    assert(fitter->fObjFunc);
00657 //    f = (*(fitter.fObjFunc) )(x);
00658 // }
00659 
00660 void TBackCompFitter::ReCreateMinimizer() { 
00661    // Recreate a minimizer instance using the function and data 
00662    // set objective function in minimizers function to re-create FCN from stored data object and fit options
00663    assert(fFitData.get());
00664 
00665    // case of standard fits (not made fia Fitter::FitFCN) 
00666    if (fFitter->Result().FittedFunction() != 0) {
00667 
00668       if (fModelFunc) delete fModelFunc; 
00669       fModelFunc =  dynamic_cast<ROOT::Math::IParamMultiFunction *>((fFitter->Result().FittedFunction())->Clone());
00670       assert(fModelFunc);
00671 
00672       // create fcn functions, should consider also gradient case
00673       const ROOT::Fit::BinData * bindata = dynamic_cast<const ROOT::Fit::BinData *>(fFitData.get()); 
00674       if (bindata) { 
00675          if (GetFitOption().Like ) 
00676             fObjFunc = new ROOT::Fit::PoissonLikelihoodFCN<ROOT::Math::IMultiGenFunction>(*bindata, *fModelFunc);
00677          else
00678             fObjFunc = new ROOT::Fit::Chi2FCN<ROOT::Math::IMultiGenFunction>(*bindata, *fModelFunc);
00679       }
00680       else { 
00681          const ROOT::Fit::UnBinData * unbindata = dynamic_cast<const ROOT::Fit::UnBinData *>(fFitData.get()); 
00682          assert(unbindata); 
00683          fObjFunc = new ROOT::Fit::LogLikelihoodFCN<ROOT::Math::IMultiGenFunction>(*unbindata, *fModelFunc);
00684       }
00685    }
00686    
00687 
00688    // recreate the minimizer
00689    fMinimizer = fFitter->Config().CreateMinimizer(); 
00690    if (fMinimizer == 0) { 
00691       Error("SetMinimizerFunction","cannot create minimizer %s",fFitter->Config().MinimizerType().c_str() );
00692    }
00693    else {
00694       if (!fObjFunc) {
00695          Error("SetMinimizerFunction","Object Function pointer is NULL");
00696       }
00697       else 
00698          fMinimizer->SetFunction(*fObjFunc);
00699    }
00700   
00701 } 
00702 
00703 
00704 
00705 void TBackCompFitter::SetFCN(void (*fcn)(Int_t &, Double_t *, Double_t &f, Double_t *, Int_t))
00706 {
00707    // override setFCN to use the Adapter to Minuit2 FCN interface
00708    //*-*-*-*-*-*-*To set the address of the minimization function*-*-*-*-*-*-*-*
00709    //*-*          ===============================================
00710    //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
00711    fFCN = fcn;
00712    if (fObjFunc) delete fObjFunc;
00713    fObjFunc = new ROOT::Fit::FcnAdapter(fFCN);
00714    DoSetDimension(); 
00715 }
00716 
00717 // need for interactive environment
00718 
00719 
00720 // global functions needed by interpreter 
00721 
00722 
00723 //______________________________________________________________________________
00724 void InteractiveFCNm2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
00725 {
00726    //*-*-*-*-*-*-*Static function called when SetFCN is called in interactive mode
00727    //*-*          ===============================================
00728    
00729    // get method call from static instance
00730    TMethodCall *m  = (TVirtualFitter::GetFitter())->GetMethodCall();
00731    if (!m) return;
00732    
00733    Long_t args[5];
00734    args[0] = (Long_t)&npar;
00735    args[1] = (Long_t)gin;
00736    args[2] = (Long_t)&f;
00737    args[3] = (Long_t)u;
00738    args[4] = (Long_t)flag;
00739    m->SetParamPtrs(args);
00740    Double_t result;
00741    m->Execute(result);
00742 }
00743 
00744 
00745 //______________________________________________________________________________
00746 void TBackCompFitter::SetFCN(void *fcn)
00747 {
00748    //*-*-*-*-*-*-*To set the address of the minimization function*-*-*-*-*-*-*-*
00749    //*-*          ===============================================
00750    //     this function is called by CINT instead of the function above
00751    //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
00752    
00753    if (!fcn) return;
00754    
00755    const char *funcname = gCint->Getp2f2funcname(fcn);
00756    if (funcname) {
00757       fMethodCall = new TMethodCall();
00758       fMethodCall->InitWithPrototype(funcname,"Int_t&,Double_t*,Double_t&,Double_t*,Int_t");
00759    }
00760    fFCN = InteractiveFCNm2;
00761    // set the static instance (required by InteractiveFCNm)
00762    TVirtualFitter::SetFitter(this); 
00763    
00764    if (fObjFunc) delete fObjFunc;
00765    fObjFunc = new ROOT::Fit::FcnAdapter(fFCN);
00766    DoSetDimension(); 
00767 }
00768 
00769 void TBackCompFitter::SetObjFunction(ROOT::Math::IMultiGenFunction   * fcn) { 
00770    // set the objective function for fitting
00771    // Needed if fitting directly using TBackCompFitter class
00772    // The class clones a copy of the function and manages it
00773    if (fObjFunc) delete fObjFunc;
00774    fObjFunc = fcn->Clone(); 
00775 }
00776 
00777 
00778 void TBackCompFitter::DoSetDimension() { 
00779    // Private method to set dimension in objective function
00780    if (!fObjFunc) return; 
00781    ROOT::Fit::FcnAdapter * fobj = dynamic_cast<ROOT::Fit::FcnAdapter*>(fObjFunc); 
00782    assert(fobj != 0); 
00783    int ndim = fFitter->Config().ParamsSettings().size(); 
00784    if (ndim != 0) fobj->SetDimension(ndim); 
00785 }
00786 
00787 ROOT::Math::IMultiGenFunction * TBackCompFitter::GetObjFunction( ) const { 
00788    // return a pointer to the objective function (FCN) 
00789    // If fitting directly using TBackCompFitter the pointer is managed by the class,
00790    // which has been set previously when calling SetObjFunction or SetFCN
00791    // Otherwise if the class is used in the backward compatible mode (e.g. after having fitted a TH1) 
00792    // the return pointer will be valid after fitting and as long a new fit will not be done. 
00793    if (fObjFunc) return fObjFunc;
00794    return fFitter->GetFCN(); 
00795 }
00796 
00797 ROOT::Math::Minimizer * TBackCompFitter::GetMinimizer( ) const { 
00798    // return a pointer to the minimizer.  
00799    // the return pointer will be valid after fitting and as long a new fit will not be done. 
00800    // For keeping a minimizer pointer the method ReCreateMinimizer() could eventually be used  
00801    if (fMinimizer) return fMinimizer;
00802    return fFitter->GetMinimizer();
00803 }
00804 
00805 TFitResult * TBackCompFitter::GetTFitResult( ) const {
00806    // return a new copy of the TFitResult object which needs to be deleted later by the user
00807    if (!fFitter.get() ) return 0; 
00808    return new TFitResult( fFitter->Result() );
00809 }
00810 
00811 //________________________________________________________________________________
00812 bool TBackCompFitter::Scan(unsigned int ipar, TGraph * gr, double xmin, double xmax )
00813 {
00814    //     scan parameter ipar between value of xmin and xmax
00815    //     a graph must be given which will be on return filled with the scan resul 
00816    //     If the graph size is zero, a default size n = 40 will be used
00817    //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
00818 
00819 
00820    if (!gr) return false; 
00821    ROOT::Math::Minimizer * minimizer = fFitter->GetMinimizer(); 
00822    if (!minimizer) {
00823       Error("Scan","Minimizer is not available - cannot scan before fitting");
00824       return false;
00825    }
00826 
00827 
00828    unsigned int npoints = gr->GetN(); 
00829    if (npoints == 0)  { 
00830       npoints = 40; 
00831       gr->Set(npoints);
00832    }
00833    bool ret = minimizer->Scan( ipar, npoints, gr->GetX(), gr->GetY(), xmin, xmax); 
00834    if ((int) npoints < gr->GetN() ) gr->Set(npoints); 
00835    return ret; 
00836 }
00837    
00838 // bool  TBackCompFitter::Scan2D(unsigned int ipar, unsigned int jpar, TGraph2D * gr, 
00839 //                       double xmin = 0, double xmax = 0, double ymin = 0, double ymax = 0) { 
00840 //    //     scan the parameters ipar between values of [xmin,xmax] and 
00841 //    //     jpar between values of [ymin,ymax] and 
00842 //    //     a graph2D must be given which will be on return filled with the scan resul 
00843 //    //     If the graph size is zero, a default size n = 20x20 will be used
00844 //    //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
00845 
00846 //    if (!gr) return false; 
00847 //    if (!fMinimizer) {
00848 //       Error("Scan","Minimizer is not available - cannot scan before fitting");
00849 //       return false;
00850 //    }
00851 //    unsigned int npoints = gr->GetN(); 
00852 //    if (npoints == 0)  { 
00853 //       npoints = 40; 
00854 //       gr->Set(npoints);
00855 //    }
00856 //    // to be implemented 
00857 //    for (unsigned int ix = 0; ix < npoints; ++ix) {       
00858 //       return fMinimizer->Scan( ipar, npoints, gr->GetX(), gr->GetY(), xmin, xmax); 
00859 
00860 // }
00861 
00862 bool  TBackCompFitter::Contour(unsigned int ipar, unsigned int jpar, TGraph * gr, double confLevel) { 
00863    //  create a 2D contour around the minimum for the parameter ipar and jpar
00864    // if a minimum does not exist or is invalid it will return false
00865    // on exit a TGraph is filled with the contour points 
00866    // the number of contur points is determined by the size of the TGraph. 
00867    // if the size is zero a default number of points = 20 is used 
00868    // pass optionally the confidence level, default is 0.683
00869    // it is assumed that ErrorDef() defines the right error definition 
00870    // (i.e 1 sigma error for one parameter). If not the confidence level are scaled to new level
00871 
00872    if (!gr) return false; 
00873    ROOT::Math::Minimizer * minimizer = fFitter->GetMinimizer(); 
00874    if (!minimizer) {
00875       Error("Scan","Minimizer is not available - cannot scan before fitting");
00876       return false;
00877    }
00878 
00879    // get error level used for fitting
00880    double upScale = fFitter->Config().MinimizerOptions().ErrorDef();
00881 
00882    double upVal = TMath::ChisquareQuantile( confLevel, 2);  // 2 is number of parameter we do the contour
00883    
00884    // set required error definition in minimizer
00885    minimizer->SetErrorDef (upScale * upVal);    
00886 
00887    unsigned int npoints = gr->GetN(); 
00888    if (npoints == 0)  { 
00889       npoints = 40; 
00890       gr->Set(npoints);
00891    }
00892    bool ret =  minimizer->Contour( ipar, jpar, npoints, gr->GetX(), gr->GetY()); 
00893    if ((int) npoints < gr->GetN() ) gr->Set(npoints); 
00894 
00895    // restore the error level used for fitting
00896    minimizer->SetErrorDef ( upScale);
00897 
00898    return ret; 
00899 }
00900 
00901 

Generated on Tue Jul 5 14:24:15 2011 for ROOT_528-00b_version by  doxygen 1.5.1