Skip to content

Commit

Permalink
Re #10532 fixed simple tests for MARI
Browse files Browse the repository at this point in the history 
and modified ReductionWrapper according to Marcus suggestions (minor)
  • Loading branch information
@abuts
abuts committed Dec 3, 2014
1 parent 1ee8235 commit 0b13082
Show file tree
Hide file tree
Showing 3 changed files with 152 additions and 757 deletions.
266 changes: 136 additions & 130 deletions Code/Mantid/scripts/Inelastic/DirectEnergyConversion.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -565,129 +565,6 @@ def _do_mono(self, data_ws, monitor_ws, result_name, ei_guess,
result_ws *= prop_man.scale_factor
return result_ws

def get_abs_normalization_factor(self,deltaE_wkspaceName,ei_monovan):
"""get absolute normalization factor for monochromatic vanadium
Inputs:
@param: deltaE_wkspace -- the name (string) of monovan workspace, converted to energy
@param: min -- the string representing the minimal energy to integrate the spectra
@param: max -- the string representing the maximal energy to integrate the spectra
@returns the value of monovan absolute normalization factor.
deletes monovan workspace (deltaE_wkspace) if abs norm factor was calculated successfully
Detailed explanation:
The algorithm takes the monochromatic vanadium workspace normalized by WB vanadium and calculates
average modified monochromatic vanadium (MV) integral considering that the efficiency of detectors
are different and accounted for by dividing each MV value by corresponding WBV value,
the signal on a detector has poison distribution and the error for a detector is the square
root of correspondent signal on a detector.
Error for WBV considered negligibly small wrt. the error on MV
"""

propman = self.prop_man;
van_mass=propman.van_mass;
minmax = propman.monovan_integr_range;


data_ws=Integration(InputWorkspace=deltaE_wkspaceName,OutputWorkspace='van_int',RangeLower=minmax[0],RangeUpper=minmax[1],IncludePartialBins='1')

nhist = data_ws.getNumberHistograms()
# extract wb integrals for combined spectra
signal=[];
error =[];
izerc=0;
for i in range(nhist):
try:
det = data_ws.getDetector(i)
except Exception:
continue
if det.isMasked():
continue
sig = data_ws.readY(i)[0]
err = data_ws.readE(i)[0]
if sig != sig: #ignore NaN (hopefully it will mean mask some day)
continue
if ((err<=0) or (sig<=0)): # count Inf and negative||zero readings. Presence of this indicates that something went wrong
izerc+=1;
continue

signal.append(sig)
error.append(err)
#---------------- Loop finished

norm_factor={};
# Various prior probabilities.
#-------------------------------------------------------------------------
# Guess which minimizes the value sum(n_i-n)^2/Sigma_i -- this what Libisis had
signal_sum = sum(map(lambda s,e: s/e,signal,error));
weight_sum = sum(map(lambda e: 1./e, error));
norm_factor['LibISIS']=signal_sum/weight_sum;
#-------------------------------------------------------------------------
# Guess which minimizes the value sum(n_i-n)^2/Sigma_i^2
signal_sum = sum(map(lambda s,e: s/(e*e),signal,error));
weight_sum = sum(map(lambda e: 1./(e*e), error));
norm_factor['SigSq']=signal_sum/weight_sum;
#-------------------------------------------------------------------------
# Guess which assumes Poisson distribution with Err=Sqrt(signal) and calculates
# the function: N_avrg = 1/(DetEfficiency_avrg^-1)*sum(n_i*DetEfficiency_i^-1)
# where the DetEfficiency = WB_signal_i/WB_average WB_signal_i is the White Beam Vanadium
# signal on i-th detector and the WB_average -- average WB vanadium signal.
# n_i is the modified signal
signal_sum = sum(map(lambda e: e*e,error));
weight_sum = sum(map(lambda s,e: e*e/s,signal,error));
if( weight_sum==0.0):
prop_man.log("WB integral has been calculated incorrectly, look at van_int workspace: {0}".format(deltaE_wkspaceName),'error')
raise ArithmeticError("Division by 0 weight when calculating WB integrals from workspace {0}".format(deltaE_wkspaceName));
norm_factor['Poisson']=signal_sum/weight_sum;
#-------------------------------------------------------------------------
# Guess which estimates value sum(n_i^2/Sigma_i^2)/sum(n_i/Sigma_i^2) TGP suggestion from 12-2012
signal_sum = sum(map(lambda s,e: s*s/(e*e),signal,error));
weight_sum = sum(map(lambda s,e: s/(e*e),signal,error));
if( weight_sum==0.0):
prop_man.log("WB integral has been calculated incorrectly, look at van_int workspace: {0}".format(deltaE_wkspaceName),'error')
raise ArithmeticError("Division by 0 weight when calculating WB integrals from workspace {0}".format(deltaE_wkspaceName));
norm_factor['TGP']=signal_sum/weight_sum;
#
#
#integral_monovan=signal_sum /(wbVan_sum)
van_multiplier = (float(propman.van_rmm)/float(van_mass))
if ei_monovan >= 210.0:
xsection = 421 # vanadium cross-section in mBarn/sR (402 mBarn/Sr) (!!!modified to fit high energy limit?!!!)
else: # old textbook cross-section for vanadium for ei=20mEv
xsection = 400 + (ei_monovan/10)
sample_multiplier = (float(propman.sample_mass)/float(propman.sample_rmm))

scale_factor = van_multiplier*sample_multiplier/xsection;

for type,val in norm_factor.iteritems():
norm_factor[type] = val*scale_factor;

# check for NaN
if (norm_factor['LibISIS']!=norm_factor['LibISIS'])|(izerc!=0): # It is an error, print diagnostics:
if (norm_factor['LibISIS'] !=norm_factor['LibISIS']):
log_value = '\n--------> Absolute normalization factor is NaN <----------------------------------------------\n'
else:
log_value ='\n--------> Warning, Monovanadium has zero spectra <--------------------------------------------\n'
log1_value = \
"""--------> Processing workspace: {0}
--------> Monovan Integration range : min={1}, max={2} (meV)
--------> Summed: {3} spectra with total signal: {4} and error: {5}
--------> Dropped: {6} zero spectra
--------> Using mBarn/sR*fu normalization factor = {7} resulting in:
--------> Abs norm factors: LibISIS: {8}
--------> Abs norm factors: Sigma^2: {9}
--------> Abs norm factors: Poisson: {10}
--------> Abs norm factors: TGP : {11}
""".format(deltaE_wkspaceName,minmax[0],minmax[1],nhist,sum(signal),sum(error),izerc,scale_factor,
norm_factor['LibISIS'],norm_factor['SigSq'],norm_factor['Poisson'],norm_factor['TGP'])
log_value = log_value+log1_value;
propman.log(log_value,'error');
else:
DeleteWorkspace(Workspace=deltaE_wkspaceName)
DeleteWorkspace(Workspace=data_ws)
return (norm_factor['LibISIS'],norm_factor['SigSq'],norm_factor['Poisson'],norm_factor['TGP'])
#-------------------------------------------------------------------------------
def convert_to_energy_transfer(self,wb_run=None,sample_run=None,ei_guess=None,rebin=None,map_file=None,monovan_run=None,wb_for_monovan_run=None,**kwargs):
""" One step conversion of run into workspace containing information about energy transfer
Expand Down Expand Up @@ -762,7 +639,6 @@ def convert_to_energy_transfer(self,wb_run=None,sample_run=None,ei_guess=None,re
# diag the sample and detector vanadium. It will deal with hard mask only if it is set that way
if not masks_done:
prop_man.log("########### Run diagnose for sample run ##############################",'notice');
prop_man.log("########### Run diagnose for monochromatic vanadium run ##############",'notice');
masking = self.diagnose(prop_man.wb_run,prop_man.mask_run,
second_white=None,print_results=True)
header = "*** Diagnostics processed workspace with {0:d} spectra and masked {1:d} bad spectra"
Expand Down Expand Up @@ -799,7 +675,7 @@ def convert_to_energy_transfer(self,wb_run=None,sample_run=None,ei_guess=None,re
prop_man.log(header.format(nSpectra,nMaskedSpectra),'notice');
#Run the conversion first on the sample
deltaE_wkspace_sample = self.convert_to_energy(prop_man.sample_run,prop_man.incident_energy,prop_man.wb_run)


# calculate absolute units integral and apply it to the workspace
if prop_man.monovan_run != None or prop_man.mono_correction_factor != None :
Expand Down Expand Up @@ -843,14 +719,16 @@ def apply_absolute_normalization(self,deltaE_wkspace_sample,monovan_run=None,ei_
"""
# Old interface support
prop_man = self.prop_man;
prop_man.log('****** Run absolute units corrections ****************************************************','notice')

if prop_man.mono_correction_factor:
absnorm_factor=float(prop_man.mono_correction_factor)
prop_man.log('##### Using supplied workspace correction factor ######','notice')
prop_man.log('*** Using supplied workspace correction factor ******','notice')
abs_norm_factor_is='provided'
else:
mvir=prop_man.monovan_integr_range;
prop_man.log('##### Evaluate the integral from the monovan run and calculate the correction factor ######','notice')
prop_man.log(' Using absolute units vanadium integration range : [{0:8f}:{1:8f}] ######'.format(mvir[0],mvir[1]),'notice')
prop_man.log('*** Evaluating the integral from the monovan run and calculate the correction factor ******','notice')
prop_man.log(' Using absolute units vanadium integration range : [{0:8f}:{1:8f}] ******'.format(mvir[0],mvir[1]),'notice')
abs_norm_factor_is = 'calculated'
#
result_ws_name = common.create_resultname(monovan_run,prop_man.instr_name)
Expand All @@ -866,7 +744,7 @@ def apply_absolute_normalization(self,deltaE_wkspace_sample,monovan_run=None,ei_
self.prop_man.map_file = map_file

ei_monovan = deltaE_wkspace_monovan.getRun().getLogData("Ei").value
prop_man.log(' Incident energy found for monovanadium run: '+str(ei_monovan)+' meV','notice')
prop_man.log(' Incident energy found for monovanadium run: '+str(ei_monovan)+' meV','notice')


(anf_LibISIS,anf_SS2,anf_Puas,anf_TGP) = self.get_abs_normalization_factor(deltaE_wkspace_monovan.getName(),ei_monovan)
Expand All @@ -875,11 +753,139 @@ def apply_absolute_normalization(self,deltaE_wkspace_sample,monovan_run=None,ei_
.format(anf_LibISIS,anf_SS2,anf_Puas,anf_TGP),'notice')
absnorm_factor = anf_TGP;
#end
prop_man.log('*** Using {0} value : {1} of absolute units correction factor'.format(abs_norm_factor_is,absnorm_factor),'notice')
prop_man.log('*** Using {0} value : {1} of absolute units correction factor (TGP)'.format(abs_norm_factor_is,absnorm_factor),'notice')
prop_man.log('*******************************************************************************************','notice')
#absnorm_factor = 0.0245159026452
# fix old system tests with would fail otherwise in 10^-6 limits
absnorm_factor = absnorm_factor*0.024519711695583177/0.0245159026452
deltaE_wkspace_sample = deltaE_wkspace_sample/absnorm_factor;


return deltaE_wkspace_sample
def get_abs_normalization_factor(self,deltaE_wkspaceName,ei_monovan):
"""get absolute normalization factor for monochromatic vanadium
Inputs:
@param: deltaE_wkspace -- the name (string) of monovan workspace, converted to energy
@param: min -- the string representing the minimal energy to integrate the spectra
@param: max -- the string representing the maximal energy to integrate the spectra
@returns the value of monovan absolute normalization factor.
deletes monovan workspace (deltaE_wkspace) if abs norm factor was calculated successfully
Detailed explanation:
The algorithm takes the monochromatic vanadium workspace normalized by WB vanadium and calculates
average modified monochromatic vanadium (MV) integral considering that the efficiency of detectors
are different and accounted for by dividing each MV value by corresponding WBV value,
the signal on a detector has poison distribution and the error for a detector is the square
root of correspondent signal on a detector.
Error for WBV considered negligibly small wrt. the error on MV
"""

propman = self.prop_man;
van_mass=propman.van_mass;
minmax = propman.monovan_integr_range;


data_ws=Integration(InputWorkspace=deltaE_wkspaceName,OutputWorkspace='van_int',RangeLower=minmax[0],RangeUpper=minmax[1],IncludePartialBins='1')

nhist = data_ws.getNumberHistograms()
# extract wb integrals for combined spectra
signal=[];
error =[];
izerc=0;
for i in range(nhist):
try:
det = data_ws.getDetector(i)
except Exception:
continue
if det.isMasked():
continue
sig = data_ws.readY(i)[0]
err = data_ws.readE(i)[0]
if sig != sig: #ignore NaN (hopefully it will mean mask some day)
continue
if ((err<=0) or (sig<=0)): # count Inf and negative||zero readings. Presence of this indicates that something went wrong
izerc+=1;
continue

signal.append(sig)
error.append(err)
#---------------- Loop finished

norm_factor={};
# Various prior probabilities.
#-------------------------------------------------------------------------
# Guess which minimizes the value sum(n_i-n)^2/Sigma_i -- this what Libisis had
signal_sum = sum(map(lambda s,e: s/e,signal,error));
weight_sum = sum(map(lambda e: 1./e, error));
norm_factor['LibISIS']=signal_sum/weight_sum;
#-------------------------------------------------------------------------
# Guess which minimizes the value sum(n_i-n)^2/Sigma_i^2
signal_sum = sum(map(lambda s,e: s/(e*e),signal,error));
weight_sum = sum(map(lambda e: 1./(e*e), error));
norm_factor['SigSq']=signal_sum/weight_sum;
#-------------------------------------------------------------------------
# Guess which assumes Poisson distribution with Err=Sqrt(signal) and calculates
# the function: N_avrg = 1/(DetEfficiency_avrg^-1)*sum(n_i*DetEfficiency_i^-1)
# where the DetEfficiency = WB_signal_i/WB_average WB_signal_i is the White Beam Vanadium
# signal on i-th detector and the WB_average -- average WB vanadium signal.
# n_i is the modified signal
signal_sum = sum(map(lambda e: e*e,error));
weight_sum = sum(map(lambda s,e: e*e/s,signal,error));
if( weight_sum==0.0):
prop_man.log("WB integral has been calculated incorrectly, look at van_int workspace: {0}".format(deltaE_wkspaceName),'error')
raise ArithmeticError("Division by 0 weight when calculating WB integrals from workspace {0}".format(deltaE_wkspaceName));
norm_factor['Poisson']=signal_sum/weight_sum;
#-------------------------------------------------------------------------
# Guess which estimates value sum(n_i^2/Sigma_i^2)/sum(n_i/Sigma_i^2) TGP suggestion from 12-2012
signal_sum = sum(map(lambda s,e: s*s/(e*e),signal,error));
weight_sum = sum(map(lambda s,e: s/(e*e),signal,error));
if( weight_sum==0.0):
prop_man.log("WB integral has been calculated incorrectly, look at van_int workspace: {0}".format(deltaE_wkspaceName),'error')
raise ArithmeticError("Division by 0 weight when calculating WB integrals from workspace {0}".format(deltaE_wkspaceName));
norm_factor['TGP']=signal_sum/weight_sum;
#
#
#integral_monovan=signal_sum /(wbVan_sum)
van_multiplier = (float(propman.van_rmm)/float(van_mass))
if ei_monovan >= 210.0:
xsection = 421 # vanadium cross-section in mBarn/sR (402 mBarn/Sr) (!!!modified to fit high energy limit?!!!)
else: # old textbook cross-section for vanadium for ei=20mEv
xsection = 400 + (ei_monovan/10)
sample_multiplier = (float(propman.sample_mass)/float(propman.sample_rmm))

scale_factor = van_multiplier*sample_multiplier/xsection;

for type,val in norm_factor.iteritems():
norm_factor[type] = val*scale_factor;

# check for NaN
if (norm_factor['LibISIS']!=norm_factor['LibISIS'])|(izerc!=0): # It is an error, print diagnostics:
if (norm_factor['LibISIS'] !=norm_factor['LibISIS']):
log_value = '\n--------> Absolute normalization factor is NaN <----------------------------------------------\n'
else:
log_value ='\n--------> Warning, Monovanadium has zero spectra <--------------------------------------------\n'
log1_value = \
"""--------> Processing workspace: {0}
--------> Monovan Integration range : min={1}, max={2} (meV)
--------> Summed: {3} spectra with total signal: {4} and error: {5}
--------> Dropped: {6} zero spectra
--------> Using mBarn/sR*fu normalization factor = {7} resulting in:
--------> Abs norm factors: LibISIS: {8}
--------> Abs norm factors: Sigma^2: {9}
--------> Abs norm factors: Poisson: {10}
--------> Abs norm factors: TGP : {11}
""".format(deltaE_wkspaceName,minmax[0],minmax[1],nhist,sum(signal),sum(error),izerc,scale_factor,
norm_factor['LibISIS'],norm_factor['SigSq'],norm_factor['Poisson'],norm_factor['TGP'])
log_value = log_value+log1_value;
propman.log(log_value,'error');
else:
DeleteWorkspace(Workspace=deltaE_wkspaceName)
DeleteWorkspace(Workspace=data_ws)
return (norm_factor['LibISIS'],norm_factor['SigSq'],norm_factor['Poisson'],norm_factor['TGP'])
#-------------------------------------------------------------------------------


def convert_to_energy(self, mono_run, ei, white_run=None, mono_van=None,\
abs_ei=None, abs_white_run=None, save_path=None, Tzero=None, \
Expand Down

0 comments on commit 0b13082

Please sign in to comment.