Source code for fudge.gnds.productData.distributions.uncorrelated
# <<BEGIN-copyright>>
# <<END-copyright>>
"""Uncorrelated double differential distribution classes."""
from __future__ import print_function
import math
from fudge.core.math import fudgemath
from PoPs import IDs as IDsPoPsModule
import xData.standards as standardsModule
import xData.ancestry as ancestryModule
import angular as angularModule
import energy as energyModule
from . import base as baseModule
__metaclass__ = type
[docs]class subform( ancestryModule.ancestry ) :
ancestryMembers = ( 'data', )
def __init__( self, data, dataSubform ) :
if( not isinstance( data, dataSubform ) ) : raise TypeError( "Needed %s distribution subform" % self.moniker )
ancestryModule.ancestry.__init__( self )
self.data = data
self.data.setAncestor( self )
@property
def productFrame( self ) :
return( self.ancestor.productFrame )
[docs] def convertUnits( self, unitMap ) :
"See documentation for reactionSuite.convertUnits."
self.data.convertUnits( unitMap )
[docs] def toXMLList( self, indent = '', **kwargs ) :
indent2 = indent + kwargs.get( 'incrementalIndent', ' ' )
xmlStringList = [ "%s<%s>" % ( indent, self.moniker ) ]
xmlStringList += self.data.toXMLList( indent2, **kwargs )
xmlStringList[-1] += "</%s>" % self.moniker
return( xmlStringList )
[docs] def to_xs_pdf_cdf1d( self, style, tempInfo, indent ) :
data = self.data.to_xs_pdf_cdf1d( style, tempInfo, indent )
if( data is None ) : return( None )
return( self.__class__( data ) )
[docs]class angularSubform( subform ) :
moniker = 'angular'
def __init__( self, data ) :
subform.__init__( self, data, angularModule.subform )
[docs]class energySubform( subform ) :
moniker = 'energy'
def __init__( self, data ) :
subform.__init__( self, data, energyModule.subform )
[docs]class form( baseModule.form ) :
"""
Contains uncorrelated distributions for outgoing angle and outgoing energy. Use when the correlations
between angle and energy are unknown.
"""
moniker = 'uncorrelated'
subformAttributes = ( 'angularSubform', 'energySubform' )
ancestryMembers = subformAttributes
def __init__( self, label, productFrame, _angularSubform, _energySubform ) :
if( not isinstance( _angularSubform, angularSubform ) ) : raise TypeError( "Needed angular distribution subform, got %s" %str(_angularSubform.__class__) )
if( not isinstance( _energySubform, energySubform ) ) : raise TypeError( "Needed energy distribution subform, got %s" % str(_energySubform.__class__) )
baseModule.form.__init__( self, label, productFrame, ( _angularSubform, _energySubform ) )
[docs] def calculateAverageProductData( self, style, indent = '', **kwargs ) :
if( isinstance( self.angularSubform.data, angularModule.regions2d ) ) :
raise Exception( 'regions2d angular is currently not supported' )
if( isinstance( self.energySubform.data, energyModule.regions2d ) ) :
aveEnergies = []
aveMomenta = []
EMax = kwargs['EMax']
for i1, energySubform in enumerate( self.energySubform.data ) :
kwargs['EMax'] = energySubform.domainMax
if( i1 == ( len( self.energySubform.data ) - 1 ) ) : kwargs['EMax'] = EMax
aveEnergy, aveMomentum = calculateAverageProductData( self.productFrame, self.angularSubform.data, energySubform,
style, indent, **kwargs )
aveEnergies.append( aveEnergy )
aveMomenta.append( aveMomentum )
kwargs['EMin'] = kwargs['EMax']
return( aveEnergies, aveMomenta )
aveEnergy, aveMomentum = calculateAverageProductData( self.productFrame, self.angularSubform.data, self.energySubform.data, style, indent, **kwargs )
return( [ aveEnergy ], [ aveMomentum ] )
[docs] def getSpectrumAtEnergy( self, energy ) :
"""Returns the energy spectrum for self at projectile energy."""
return( self.energySubform.data.getSpectrumAtEnergy( energy ) )
[docs] def findEntity( self, entityName, attribute = None, value = None ):
"""
Overrides ancestry.findEntity. uncorrelated contains both energy and angular form,
need to check both to find desired entity
"""
if self.energySubform.moniker == entityName: return self.energySubform
elif self.angularSubform.moniker == entityName: return self.angularSubform
return ancestryModule.ancestry.findEntity( self, entityName, attribute, value )
[docs] def processMC_cdf( self, style, tempInfo, indent ) :
_angularSubform = self.angularSubform.to_xs_pdf_cdf1d( style, tempInfo, indent )
_energySubform = self.energySubform.to_xs_pdf_cdf1d( style, tempInfo, indent )
if( ( _angularSubform is not None ) or ( _energySubform is not None ) ) :
if( _angularSubform is None ) : _angularSubform = self.angularSubform.copy( )
if( _energySubform is None ) : _energySubform = self.energySubform.copy( )
_form = form( style.label, self.productFrame, _angularSubform, _energySubform )
else :
_form = None
return( _form )
[docs] def processMultiGroup( self, style, tempInfo, indent ) :
from fudge.processing import group as groupModule
from fudge.processing.deterministic import transferMatrices as transferMatricesModule
verbosity = tempInfo['verbosity']
indent2 = indent + tempInfo['incrementalIndent']
productLabel = tempInfo['productLabel']
angularSubform = self.angularSubform.data
energySubform = self.energySubform.data
energyUnit = tempInfo['incidentEnergyUnit']
massUnit = energyUnit + '/c**2'
if( verbosity > 2 ) : print('%sGrouping %s' % (indent, self.moniker))
crossSection = tempInfo['crossSection']
product = tempInfo['product']
if( isinstance( energySubform, energyModule.discreteGamma ) ) :
if( product.id == IDsPoPsModule.photon ) :
Ep = float( energySubform.value )
TM_1, TM_E = transferMatricesModule.discreteGammaAngularData( style, tempInfo, Ep, crossSection,
angularSubform, multiplicity = product.multiplicity.evaluated,
comment = tempInfo['transferMatrixComment'] + ' outgoing data for %s' % productLabel )
else :
raise Exception( 'See Bret' )
else :
if( isinstance( energySubform, energyModule.NBodyPhaseSpace ) ) :
totalMass = energySubform.numberOfProductsMasses.getValueAs( massUnit )
Q = tempInfo['reaction'].getQ( energyUnit, final = False )
TM_1, TM_E = transferMatricesModule.NBodyPhaseSpace( style, tempInfo, crossSection,
energySubform.numberOfProducts, totalMass, Q, tempInfo['multiplicity'],
comment = tempInfo['transferMatrixComment'] + ' outgoing data for %s' % productLabel )
else :
TM_1, TM_E = transferMatricesModule.uncorrelated_EMuP_EEpP_TransferMatrix( style, tempInfo, crossSection,
self.productFrame, angularSubform, energySubform, tempInfo['multiplicity'],
comment = tempInfo['transferMatrixComment'] + ' outgoing data for %s' % productLabel )
return( groupModule.TMs2Form( style, tempInfo, TM_1, TM_E ) )
[docs] def toPointwise_withLinearXYs( self, **kwargs ) :
import angularEnergy
angularSubform = self.angularSubform.toPointwise_withLinearXYs( **kwargs )
energySubform = self.energySubform.toPointwise_withLinearXYs( **kwargs )
grid = angularSubform.domainGrid
for e in energySubform.domainGrid :
if( e not in grid ) : grid.append( e )
grid.sort( )
axes = angularEnergy.defaultAxes( energyUnit = energySubform.axes[-1].unit )
f_E_mu_Ep = angularEnergy.XYs3d( axes = axes )
for e in grid :
f_mu_Ep = angularEnergy.XYs2d( value = e )
af = angularSubform.getAtEnergy( e ) # FIXME why doesn't getSpectrumAtEnergy work here?
ef = energySubform.getSpectrumAtEnergy( e )
for mu, P in af :
efp = P * ef
efp.__class__ = angularEnergy.XYs1d # FIXME better way to assign class?
efp.value = mu
f_mu_Ep.append( efp )
f_E_mu_Ep.append( f_mu_Ep )
return( f_E_mu_Ep )
[docs] def toXMLList( self, indent = '', **kwargs ) :
indent2 = indent + kwargs.get( 'incrementalIndent', ' ' )
attributeStr = ''
if( self.label is not None ) : attributeStr += ' label="%s"' % self.label
if( self.productFrame is not None ) : attributeStr += ' productFrame="%s"' % self.productFrame
xmlString = [ '%s<%s%s>' % ( indent, self.moniker, attributeStr ) ]
xmlString += self.angularSubform.toXMLList( indent2, **kwargs )
xmlString += self.energySubform.toXMLList( indent2, **kwargs )
xmlString[-1] += '</%s>' % self.moniker
return( xmlString )
[docs] @staticmethod
def parseXMLNode( element, xPath, linkData ) :
"""Translate <uncorrelated> element from xml. Must contain one <angular> and one <energy> component."""
xPath.append( element.tag )
for child in element :
if( child.tag == angularModule.form.moniker ) :
_angularSubform = angularSubform( angularModule.form.parseXMLNode( child[0], xPath, linkData ) )
elif( child.tag == energyModule.form.moniker ) :
_energySubform = energySubform( energyModule.form.parseXMLNode( child[0], xPath, linkData ) )
else :
raise ValueError( 'unsupported tag = "%s"' % child.tag )
uncorrelated = form( element.get( 'label' ), element.get( 'productFrame' ), _angularSubform, _energySubform )
xPath.pop( )
return( uncorrelated )
[docs]def calculateAverageProductData( productFrame, angularSubform, energySubform, style, indent, **kwargs ) :
def fixLimits( EMin, EMax, Es ) :
if( EMin not in Es ) : Es.append( EMin )
if( EMax not in Es ) : Es.append( EMax )
Es.sort( )
while( Es[0] < EMin ) : del Es[0]
while( Es[-1] > EMax ) : del Es[-1]
def calculateAverageEnergy( self, Ein ) :
averageEp = self.energySubform.averageEp( Ein )
if( self.productFrame == standardsModule.frames.centerOfMassToken ) :
A_Ein = self.massRatio * Ein
if( not( self.angularSubform.isIsotropic( ) ) ) : averageEp += 2 * math.sqrt( A_Ein * averageEp ) * self.angularSubform.averageMu( Ein )
averageEp += A_Ein
averageEp *= self.multiplicity.evaluate( Ein )
return( averageEp )
def calculateAverageMomentum( self, Ein ) :
pp, averageMu = 0., self.angularSubform.averageMu( E )
if( averageMu != 0. ) : pp = energySubform.sqrtEp_AverageAtE( E ) * averageMu
if( self.productFrame == standardsModule.frames.centerOfMassToken ) : pp += math.sqrt( self.massRatio * Ein )
return( multiplicity.evaluate( E ) * math.sqrt( 2. * self.productMass ) * pp )
def calculateAverageMomentumForPhoton( self, Ein ) :
muAverage = angularSubform.averageMu( E )
if( muAverage == 0. ) : return( 0. )
return( multiplicity.evaluate( E ) * energySubform.averageEp( E ) * muAverage )
class calculateDepositionInfo :
def __init__( self, productFrame, productMass, massRatio, angularSubform, energySubform, multiplicity ) :
self.productFrame = productFrame
self.productMass = productMass
self.massRatio = massRatio
self.angularSubform = angularSubform
self.energySubform = energySubform
self.multiplicity = multiplicity
def evaluateAtX( self, E ) :
return( self._evaluateAtX( self, E ) )
def setEvaluateAtX( self, evaluateAtX ) :
self._evaluateAtX = evaluateAtX
def setTolerances( self, relativeTolerance, absoluteTolerance ) :
self.relativeTolerance = relativeTolerance
self.absoluteTolerance = absoluteTolerance
class NBodyPhaseSpace :
def __init__( self, energySubform, massUnit, projectileMass, targetMass, productMass, Q ) :
self.energySubform = energySubform
self.massUnit = massUnit
self.projectileMass = projectileMass
self.targetMass = targetMass
self.productMass = productMass
self.Q = Q
def averageEp( self, Ein ) :
return( self.energySubform.averageEp( Ein, self.massUnit, self.projectileMass, self.targetMass, self.productMass, self.Q ) )
def getEnergyArray( self, EMin = None, EMax = None ) :
return( [ EMin, EMax ] )
energyUnit = kwargs['incidentEnergyUnit']
massUnit = kwargs['massUnit']
multiplicity = kwargs['multiplicity'] # BRB - FIXME; Handle gamma data with 1 point for multiplicity weight until it is fixed.
energyAccuracy = kwargs['energyAccuracy']
momentumAccuracy = kwargs['momentumAccuracy']
projectileMass = kwargs['projectileMass']
targetMass = kwargs['targetMass']
product = kwargs['product']
productMass = kwargs['productMass']
EMin = kwargs['EMin']
EMax = kwargs['EMax']
massRatio = projectileMass * productMass / ( projectileMass + targetMass )**2
if( product.id == IDsPoPsModule.photon ) : productFrame = standardsModule.frames.labToken # All gamma data treated as in lab frame.
if( isinstance( energySubform, energyModule.NBodyPhaseSpace ) ) :
Q = kwargs['reaction'].getQ( energyUnit, final = False )
energySubform = NBodyPhaseSpace( energySubform, massUnit, projectileMass, targetMass, productMass, Q )
calculationData = calculateDepositionInfo( productFrame, productMass, massRatio, angularSubform, energySubform, multiplicity )
Es = energySubform.getEnergyArray( )
if( Es[0] is None ) : Es[0] = EMin
if( Es[-1] is None ) : Es[-1] = EMax
if( EMin < Es[0] ) : EMin = Es[0]
if( EMax > Es[-1] ) : EMax = Es[-1]
Es = sorted( set( energySubform.getEnergyArray( EMin, EMax ) + multiplicity.domainGrid ) )
fixLimits( EMin, EMax, Es )
calculationData.setEvaluateAtX( calculateAverageEnergy )
aveEnergy = [ [ E, calculationData.evaluateAtX( E ) ] for E in Es ]
absoluteTolerance = 1e-3 * energyAccuracy * max( [ Ep for E, Ep in aveEnergy ] )
calculationData.setTolerances( energyAccuracy, absoluteTolerance )
aveEnergy = fudgemath.thickenXYList( aveEnergy, calculationData )
if( isinstance( angularSubform, angularModule.isotropic2d ) and ( productFrame == standardsModule.frames.labToken ) ) :
aveMomentum = [ [ aveEnergy[0][0], 0. ], [ aveEnergy[-1][0], 0. ] ]
else :
if( product.id == IDsPoPsModule.photon ) :
calculationData.setEvaluateAtX( calculateAverageMomentumForPhoton )
else :
calculationData.setEvaluateAtX( calculateAverageMomentum )
Es = sorted( set( Es + angularSubform.getEnergyArray( EMin, EMax ) ) )
fixLimits( EMin, EMax, Es )
aveMomentum = [ [ E, calculationData.evaluateAtX( E ) ] for E in Es ]
absoluteTolerance = 1e-3 * momentumAccuracy * max( [ pp for E, pp in aveMomentum ] )
calculationData.setTolerances( momentumAccuracy, absoluteTolerance )
aveMomentum = fudgemath.thickenXYList( aveMomentum, calculationData )
return( aveEnergy, aveMomentum )