def profile(self, name, parName, logemin=None, logemax=None,
                reoptimize=False,
                xvals=None, npts=None, savestate=True, **kwargs):
        """Profile the likelihood for the given source and parameter.

        Parameters
        ----------

        name : str
           Source name.

        parName : str
           Parameter name.

        reoptimize : bool
           Re-fit nuisance parameters at each step in the scan.  Note
           that enabling this option will only re-fit parameters that
           were free when the method was executed.

        Returns
        -------

        lnlprofile : dict
           Dictionary containing results of likelihood scan.

        """
        # Find the source
        name = self.roi.get_source_by_name(name).name

        par = self.like.normPar(name)
        parName = self.like.normPar(name).getName()
        idx = self.like.par_index(name, parName)
        bounds = self.like.model[idx].getBounds()
        value = self.like.model[idx].getValue()
        loge_bounds = self.loge_bounds

        optimizer = kwargs.get('optimizer', self.config['optimizer'])

        if savestate:
            saved_state = self._latch_state()

        # If parameter is fixed temporarily free it
        par.setFree(True)
        if optimizer['optimizer'] == 'NEWTON':
            self._create_fitcache()

        if logemin is not None or logemax is not None:
            loge_bounds = self.set_energy_range(logemin, logemax)
        else:
            loge_bounds = self.loge_bounds

        loglike0 = -self.like()

        if xvals is None:

            err = par.error()
            val = par.getValue()
            if err <= 0 or val <= 3 * err:
                xvals = 10 ** np.linspace(-2.0, 2.0, 51)
                if val < xvals[0]:
                    xvals = np.insert(xvals, val, 0)
            else:
                xvals = np.linspace(0, 1, 25)
                xvals = np.concatenate((-1.0 * xvals[1:][::-1], xvals))
                xvals = val * 10 ** xvals

        if np.isnan(xvals).any():
            raise RuntimeError(
                "Parameter scan points for %s::%s include infinite value." % (name, parName))

        # Update parameter bounds to encompass scan range
        try:
            self.like[idx].setBounds(min(min(xvals), value, bounds[0]),
                                     max(max(xvals), value, bounds[1]))
        except RuntimeError:
            self.logger.warning(
                "Caught failure on setBounds for %s::%s." % (name, parName))

        o = {'xvals': xvals,
             'npred': np.zeros(len(xvals)),
             'npred_wt': np.zeros(len(xvals)),
             'dnde': np.zeros(len(xvals)),
             'flux': np.zeros(len(xvals)),
             'eflux': np.zeros(len(xvals)),
             'dloglike': np.zeros(len(xvals)),
             'loglike': np.zeros(len(xvals))
             }

        if reoptimize and hasattr(self.like.components[0].logLike,
                                  'setUpdateFixedWeights'):

            for c in self.components:
                c.like.logLike.setUpdateFixedWeights(False)

        for i, x in enumerate(xvals):

            try:
                self.like[idx] = x
            except RuntimeError:
                self.logger.warning(
                    "Caught failure on set for %s::%s: %.2f" % (name, parName, x))

            if self.like.nFreeParams() > 1 and reoptimize:
                # Only reoptimize if not all frozen
                self.like.freeze(idx)
                fit_output = self._fit(errors=False, **optimizer)
                loglike1 = fit_output['loglike']
                self.like.thaw(idx)
            else:
                loglike1 = -self.like()

            flux = self.like[name].flux(10 ** loge_bounds[0],
                                        10 ** loge_bounds[1])
            eflux = self.like[name].energyFlux(10 ** loge_bounds[0],
                                               10 ** loge_bounds[1])
            prefactor = self.like[idx]

            o['dloglike'][i] = loglike1 - loglike0
            o['loglike'][i] = loglike1
            o['dnde'][i] = prefactor.getTrueValue()
            o['flux'][i] = flux
            o['eflux'][i] = eflux

            cs = self.model_counts_spectrum(name,
                                            loge_bounds[0],
                                            loge_bounds[1], summed=True)
            o['npred'][i] += np.sum(cs)
            cs_wt = self.model_counts_spectrum(name,
                                               loge_bounds[0],
                                               loge_bounds[1], summed=True, weighted=True)
            o['npred_wt'][i] += np.sum(cs_wt)

        self.like[idx] = value

        if reoptimize and hasattr(self.like.components[0].logLike,
                                  'setUpdateFixedWeights'):

            for c in self.components:
                c.like.logLike.setUpdateFixedWeights(True)

        # Restore model parameters to original values
        if savestate:
            saved_state.restore()

        self.like[idx].setBounds(*bounds)
        if logemin is not None or logemax is not None:
            self.set_energy_range(*loge_bounds)

        return o