00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 #include "Minuit2/MnParameterScan.h"
00011 #include "Minuit2/FCNBase.h"
00012 
00013 namespace ROOT {
00014 
00015    namespace Minuit2 {
00016 
00017 
00018 MnParameterScan::MnParameterScan(const FCNBase& fcn, const MnUserParameters& par) : fFCN(fcn), fParameters(par), fAmin(fcn(par.Params())) {}
00019 
00020 MnParameterScan::MnParameterScan(const FCNBase& fcn, const MnUserParameters& par, double fval) : fFCN(fcn), fParameters(par), fAmin(fval) {}
00021 
00022 std::vector<std::pair<double, double> > MnParameterScan::operator()(unsigned int par, unsigned int maxsteps, double low, double high) {
00023    
00024    
00025    
00026    std::vector<std::pair<double, double> > result; result.reserve(maxsteps+1);
00027    std::vector<double> params = fParameters.Params();
00028    result.push_back(std::pair<double, double>(params[par], fAmin));
00029    
00030    if(low > high) return result;
00031    if(maxsteps < 2) return result;
00032    
00033    if(low == 0. && high == 0.) {
00034       low = params[par] - 2.*fParameters.Error(par);
00035       high = params[par] + 2.*fParameters.Error(par);
00036    }
00037    
00038    if(low == 0. && high == 0. && fParameters.Parameter(par).HasLimits()) {
00039       if(fParameters.Parameter(par).HasLowerLimit())
00040          low = fParameters.Parameter(par).LowerLimit();
00041       if(fParameters.Parameter(par).HasUpperLimit())
00042          high = fParameters.Parameter(par).UpperLimit();
00043    }
00044    
00045    if(fParameters.Parameter(par).HasLimits()) {
00046       if(fParameters.Parameter(par).HasLowerLimit())
00047          low = std::max(low, fParameters.Parameter(par).LowerLimit());
00048       if(fParameters.Parameter(par).HasUpperLimit())
00049          high = std::min(high, fParameters.Parameter(par).UpperLimit());
00050    }
00051    
00052    double x0 = low;
00053    double stp = (high - low)/double(maxsteps - 1);
00054    for(unsigned int i = 0; i < maxsteps; i++) {
00055       params[par] = x0 + double(i)*stp;
00056       double fval = fFCN(params);
00057       if(fval < fAmin) {
00058          fParameters.SetValue(par, params[par]);
00059          fAmin = fval;
00060       }
00061       result.push_back(std::pair<double, double>(params[par], fval));
00062    }
00063    
00064    return result;
00065 }
00066 
00067    }  
00068 
00069 }