diff --git a/src/mvesuvio/analysis_fitting.py b/src/mvesuvio/analysis_fitting.py index d866764..3f955e9 100644 --- a/src/mvesuvio/analysis_fitting.py +++ b/src/mvesuvio/analysis_fitting.py @@ -15,7 +15,7 @@ repoPath = Path(__file__).absolute().parent # Path to the repository -def fitInYSpaceProcedure(yFitIC, IC, wsTOF): +def fitInYSpaceProcedure(IC, wsTOF): try: wsTOF = mtd[wsTOF] @@ -23,31 +23,31 @@ def fitInYSpaceProcedure(yFitIC, IC, wsTOF): print(f"Workspace to fit {wsTOF} not found.") return - wsResSum, wsRes = calculateMantidResolutionFirstMass(IC, yFitIC, wsTOF) + wsResSum, wsRes = calculateMantidResolutionFirstMass(IC, wsTOF) wsTOFMass0 = subtract_profiles_except_lightest(IC, wsTOF) - if yFitIC.subtract_calculated_fse_from_data: + if IC.subtract_calculated_fse_from_data: wsFSEMass0 = find_ws_name_fse_first_mass(IC) wsTOFMass0 = Minus(wsTOFMass0, wsFSEMass0, OutputWorkspace=wsTOFMass0.name() + "_fse") - wsJoY, wsJoYAvg = ySpaceReduction(wsTOFMass0, IC.masses[0], yFitIC, IC) + wsJoY, wsJoYAvg = ySpaceReduction(wsTOFMass0, IC) - if yFitIC.do_symmetrisation: + if IC.do_symmetrisation: wsJoYAvg = symmetrizeWs(wsJoYAvg) - fitProfileMinuit(yFitIC, wsJoYAvg, wsResSum) - fitProfileMantidFit(yFitIC, wsJoYAvg, wsResSum) + fitProfileMinuit(IC, wsJoYAvg, wsResSum) + fitProfileMantidFit(IC, wsJoYAvg, wsResSum) printYSpaceFitResults(wsJoYAvg.name()) - yfitResults = ResultsYFitObject(IC, yFitIC, wsTOF.name(), wsJoYAvg.name()) + yfitResults = ResultsYFitObject(IC, wsTOF.name(), wsJoYAvg.name()) yfitResults.save() - if yFitIC.do_global_fit: - runGlobalFit(wsJoY, wsRes, IC, yFitIC) + if IC.do_global_fit: + runGlobalFit(wsJoY, wsRes, IC) - save_workspaces(yFitIC) + save_workspaces(IC) return yfitResults @@ -67,7 +67,7 @@ def find_ws_name_fse_first_mass(ic): return ws_names_fse[np.argmin(ws_masses)] -def calculateMantidResolutionFirstMass(IC, yFitIC, ws): +def calculateMantidResolutionFirstMass(IC, ws): mass = IC.masses[0] resName = ws.name() + "_Resolution" @@ -77,7 +77,7 @@ def calculateMantidResolutionFirstMass(IC, yFitIC, ws): ) Rebin( InputWorkspace="tmp", - Params=yFitIC.range_for_rebinning_in_y_space, + Params=IC.range_for_rebinning_in_y_space, OutputWorkspace="tmp", ) @@ -124,19 +124,20 @@ def switchFirstTwoAxis(A): return np.stack(np.split(A, len(A), axis=0), axis=2)[0] -def ySpaceReduction(wsTOF, mass0, yFitIC, ic): +def ySpaceReduction(wsTOF, ic): """Seperate procedures depending on masking specified.""" + mass0 = ic.masses[0] ws_name_lightest_mass = ic.name + '_' + str(ic.number_of_iterations_for_corrections) + '_' + str(min(ic.masses)) + '_ncp' ncp = mtd[ws_name_lightest_mass].extractY() - rebinPars = yFitIC.range_for_rebinning_in_y_space + rebinPars = ic.range_for_rebinning_in_y_space if np.any(np.all(wsTOF.extractY() == 0, axis=0)): # Masked columns present - if yFitIC.mask_zeros_with == "nan": + if ic.mask_zeros_with == "nan": # Build special workspace to store accumulated points wsJoY = convertToYSpace(wsTOF, mass0) - xp = buildXRangeFromRebinPars(yFitIC) + xp = buildXRangeFromRebinPars(ic) wsJoYB = dataXBining( wsJoY, xp ) # Unusual ws with several dataY points per each dataX point @@ -1319,7 +1320,7 @@ def printYSpaceFitResults(wsJoYName): class ResultsYFitObject: - def __init__(self, ic, yFitIC, wsFinalName, wsYSpaceAvgName): + def __init__(self, ic, wsFinalName, wsYSpaceAvgName): # Extract most relevant information from ws wsFinal = mtd[wsFinalName] wsResSum = mtd[wsFinalName + "_Resolution_Sum"] @@ -1339,19 +1340,19 @@ def __init__(self, ic, yFitIC, wsFinalName, wsYSpaceAvgName): poptList = [] perrList = [] try: - wsFitMinuit = mtd[wsJoYAvg.name() + "_minuit_" + yFitIC.fitting_model + "_Parameters"] + wsFitMinuit = mtd[wsJoYAvg.name() + "_minuit_" + ic.fitting_model + "_Parameters"] poptList.append(wsFitMinuit.column("Value")) perrList.append(wsFitMinuit.column("Error")) except: pass try: - wsFitLM = mtd[wsJoYAvg.name() + "_lm_" + yFitIC.fitting_model + "_Parameters"] + wsFitLM = mtd[wsJoYAvg.name() + "_lm_" + ic.fitting_model + "_Parameters"] poptList.append(wsFitLM.column("Value")) perrList.append(wsFitLM.column("Error")) except: pass try: - wsFitSimplex = mtd[wsJoYAvg.name() + "_simplex_" + yFitIC.fitting_model + "_Parameters"] + wsFitSimplex = mtd[wsJoYAvg.name() + "_simplex_" + ic.fitting_model + "_Parameters"] poptList.append(wsFitSimplex.column("Value")) perrList.append(wsFitSimplex.column("Error")) except: @@ -1371,7 +1372,7 @@ def __init__(self, ic, yFitIC, wsFinalName, wsYSpaceAvgName): self.perr = perr self.savePath = ic.ySpaceFitSavePath - self.fitting_model= yFitIC.fitting_model + self.fitting_model= ic.fitting_model def save(self): np.savez( @@ -1387,7 +1388,7 @@ def save(self): ) -def runGlobalFit(wsYSpace, wsRes, IC, yFitIC): +def runGlobalFit(wsYSpace, wsRes, IC): print("\nRunning GLobal Fit ...\n") dataX, dataY, dataE, dataRes, instrPars = extractData(wsYSpace, wsRes, IC) @@ -1395,15 +1396,15 @@ def runGlobalFit(wsYSpace, wsRes, IC, yFitIC): dataX, dataY, dataE, dataRes, instrPars ) - idxList = groupDetectors(instrPars, yFitIC) + idxList = groupDetectors(instrPars, IC) dataX, dataY, dataE, dataRes = avgWeightDetGroups( - dataX, dataY, dataE, dataRes, idxList, yFitIC + dataX, dataY, dataE, dataRes, idxList, IC ) - if yFitIC.do_symmetrisation: + if IC.do_symmetrisation: dataY, dataE = weightedSymArr(dataY, dataE) - model, defaultPars, sharedPars = selectModelAndPars(yFitIC.fitting_model) + model, defaultPars, sharedPars = selectModelAndPars(IC.fitting_model) totCost = 0 for i, (x, y, yerr, res) in enumerate(zip(dataX, dataY, dataE, dataRes)): @@ -1424,12 +1425,12 @@ def runGlobalFit(wsYSpace, wsRes, IC, yFitIC): for i in range(len(dataY)): # Set limits for unshared parameters m.limits["A" + str(i)] = (0, np.inf) - if yFitIC.fitting_model== "doublewell": + if IC.fitting_model== "doublewell": m.limits["d"] = (0, np.inf) # Shared parameters m.limits["R"] = (0, np.inf) t0 = time.time() - if yFitIC.fitting_model== "gauss": + if IC.fitting_model== "gauss": m.simplex() m.migrad() @@ -1477,8 +1478,8 @@ def constr(*pars): create_table_for_global_fit_parameters(wsYSpace.name(), m, chi2) - if yFitIC.show_plots: - plotGlobalFit(dataX, dataY, dataE, m, totCost, wsYSpace.name(), yFitIC) + if IC.show_plots: + plotGlobalFit(dataX, dataY, dataE, m, totCost, wsYSpace.name(), IC) # Pass into array to store values in variable return np.array(m.values), np.array(m.errors) diff --git a/src/mvesuvio/config/analysis_inputs.py b/src/mvesuvio/config/analysis_inputs.py index 0bc4039..c7605f4 100644 --- a/src/mvesuvio/config/analysis_inputs.py +++ b/src/mvesuvio/config/analysis_inputs.py @@ -11,8 +11,8 @@ class SampleParameters: @dataclass class BackwardAnalysisInputs(SampleParameters): - run_this_scattering_type = False - fit_in_y_space = False + run_this_scattering_type = True + fit_in_y_space = True runs = "43066-43076" empty_runs = "41876-41923" @@ -45,6 +45,29 @@ class BackwardAnalysisInputs(SampleParameters): do_gamma_correction = False + show_plots = True + do_symmetrisation = False + subtract_calculated_fse_from_data = True + range_for_rebinning_in_y_space = "-25, 0.5, 25" # Needs to be symetric + # Fitting model options + # 'gauss': Single Gaussian + # 'gcc4': Gram-Charlier with C4 parameter + # 'gcc6': Gram-Charlier with C6 parameter + # 'doublewell': Double Well function + # 'gauss3D': 3-Dimensional Gaussian + fitting_model = "gauss" + run_minos = True + do_global_fit = True # Performs global fit with Minuit by default + # Number of groups of detectors to perform global (simultaneous) fit on + # Either an integer less than the number of detectors + # or option 'all', which does not form groups and fits all spectra simultaneously and individualy + number_of_global_fit_groups = 4 + # Type of masking + # 'nan': Zeros in workspace being fit are ignored + # 'ncp': Zeros in workspace being fit are replaced by the fitted neutron compton profile + mask_zeros_with = "nan" + + @dataclass class ForwardAnalysisInputs(SampleParameters): run_this_scattering_type = True @@ -81,8 +104,6 @@ class ForwardAnalysisInputs(SampleParameters): do_gamma_correction = True -@dataclass -class YSpaceFitInputs: show_plots = True do_symmetrisation = False subtract_calculated_fse_from_data = True diff --git a/src/mvesuvio/main/run_routine.py b/src/mvesuvio/main/run_routine.py index 354a9c0..57fd01e 100644 --- a/src/mvesuvio/main/run_routine.py +++ b/src/mvesuvio/main/run_routine.py @@ -32,7 +32,6 @@ def setup(self): self.bckwd_ai = ai.BackwardAnalysisInputs self.fwd_ai = ai.ForwardAnalysisInputs - self.yFitIC = ai.YSpaceFitInputs # Names of workspaces to check if they exist to skip analysis self.ws_to_fit_y_space = [] @@ -63,7 +62,8 @@ def setup(self): # TODO: sort out yfit inputs figSavePath = self.experiment_path / "figures" figSavePath.mkdir(exist_ok=True) - self.yFitIC.figSavePath = figSavePath + self.fwd_ai.figSavePath = figSavePath + self.bckwd_ai.figSavePath = figSavePath def import_from_inputs(self): @@ -95,7 +95,7 @@ def run(self): def runAnalysisFitting(self): for wsName, i_cls in zip(self.ws_to_fit_y_space, self.classes_to_fit_y_space): - self.fitting_result = fitInYSpaceProcedure(self.yFitIC, i_cls, wsName) + self.fitting_result = fitInYSpaceProcedure(i_cls, wsName) return