# <<BEGIN-copyright>>
# <<END-copyright>>
"""Energy/angular double differential distribution classes."""
import math
import xData.axes as axesModule
import xData.standards as standardsModule
import xData.series1d as series1dModule
import xData.XYs as XYsModule
import xData.multiD_XYs as multiD_XYsModule
import xData.regions as regionsModule
from . import base as baseModule
__metaclass__ = type
[docs]def defaultAxes( energyUnit ) :
axes = axesModule.axes( rank = 4 )
axes[0] = axesModule.axis( 'P(energy_out,mu|energy_in)', 0, '1/' + energyUnit )
axes[1] = axesModule.axis( 'mu', 1, '' )
axes[2] = axesModule.axis( 'energy_out', 2, energyUnit )
axes[3] = axesModule.axis( 'energy_in', 3, energyUnit )
return( axes )
[docs]class XYs1d( XYsModule.XYs1d ) : # FIXME: BRB, Should this class inherit from angular.XYs1d?
[docs] def averageMu( self ) :
allowedInterpolations = [ standardsModule.interpolation.linlinToken,
standardsModule.interpolation.flatToken ]
xys = self.changeInterpolationIfNeeded( allowedInterpolations, XYsModule.defaultAccuracy )
return( xys.integrateWithWeight_x( ) )
[docs]class Legendre( series1dModule.LegendreSeries ) :
[docs] def averageMu( self ) :
return( self.getCoefficientSafely( 1 ) )
[docs] def integrate( self ) :
return( self.getCoefficientSafely( 0 ) )
[docs]class XYs2d( multiD_XYsModule.XYs2d ) :
[docs] def averageEnergy( self ) :
integralOfMu = [ [ pdfOfMuAtEp.value, float( pdfOfMuAtEp.integrate( ) ) ] for pdfOfMuAtEp in self ]
return( float( XYsModule.XYs1d( integralOfMu, interpolation = self.interpolation ).integrateWithWeight_x( ) ) )
[docs] def averageForwardMomentum( self, sqrt_2_ProductMass ) :
averageMu = [ [ pdfOfMuAtEp.value, pdfOfMuAtEp.averageMu( ) ] for pdfOfMuAtEp in self ]
return( sqrt_2_ProductMass * float( XYsModule.XYs1d( averageMu, interpolation = self.interpolation ).integrateWithWeight_sqrt_x( ) ) )
[docs] @staticmethod
def allowedSubElements( ) :
return( ( XYs1d, Legendre ) )
[docs]class XYs3d( subform, multiD_XYsModule.XYs3d ) :
def __init__( self, **kwargs ) :
multiD_XYsModule.XYs3d.__init__( self, **kwargs )
subform.__init__( self )
[docs] def calculateAverageProductData( self, style, indent = '', **kwargs ) :
verbosity = kwargs.get( 'verbosity', 0 )
indent2 = indent + kwargs.get( 'incrementalIndent', ' ' )
multiplicity = kwargs['multiplicity']
projectileMass = kwargs['projectileMass']
targetMass = kwargs['targetMass']
productMass = kwargs['productMass']
energyAccuracy = kwargs['energyAccuracy']
momentumAccuracy = kwargs['momentumAccuracy']
productFrame = kwargs['productFrame']
sqrt_2_ProductMass = math.sqrt( 2 * productMass )
#
# FIXME, still need to fill in between energy points.
#
aveEnergy = []
aveMomentum = []
if( productFrame == standardsModule.frames.labToken ) :
for pdfOfEpMuAtE in self :
energy = pdfOfEpMuAtE.value
aveEnergy.append( [ energy, multiplicity.evaluate( energy ) * pdfOfEpMuAtE.averageEnergy( ) ] )
for pdfOfEpMuAtE in self :
energy = pdfOfEpMuAtE.value
momemtum = multiplicity.evaluate( energy ) * pdfOfEpMuAtE.averageForwardMomentum( sqrt_2_ProductMass )
if( momemtum < 1e-12 ) : momemtum = 0. # This should be less arbitrary????????
aveMomentum.append( [ energy, momemtum ] )
else : # Center-of-mass.
const1 = math.sqrt( 2 * projectileMass ) / ( projectileMass + targetMass )
for pdfOfEpMuAtE in self :
energy = pdfOfEpMuAtE.value
vCOM = const1 * math.sqrt( energy )
EpCOM = pdfOfEpMuAtE.averageEnergy( )
ppCOM = pdfOfEpMuAtE.averageForwardMomentum( sqrt_2_ProductMass )
productLabEnergy = 0.5 * productMass * vCOM * vCOM + EpCOM + vCOM * ppCOM
aveEnergy.append( [ energy, multiplicity.evaluate( energy ) * productLabEnergy ] )
for pdfOfEpMuAtE in self :
energy = pdfOfEpMuAtE.value
vCOM = const1 * math.sqrt( energy )
productLabForwardMomentum = productMass * vCOM + pdfOfEpMuAtE.averageForwardMomentum( sqrt_2_ProductMass )
aveMomentum.append( [ energy, multiplicity.evaluate( energy ) * productLabForwardMomentum ] )
return( [ aveEnergy ], [ aveMomentum ] )
[docs] def check( self, info ) :
from fudge.gnds import warning
from pqu import PQU
warnings = []
axes = axesModule.axes()
axes[0] = self.axes[-2]
axes[1] = self.axes[-1]
for index, energy_in in enumerate(self):
integral = float( XYsModule.XYs1d( [ ( eout.value, eout.coefficients[0] ) for eout in energy_in ], axes = axes ).integrate( ) )
if abs(integral-1.0) > info['normTolerance']:
warnings.append( warning.unnormalizedDistribution( PQU.PQU(energy_in.value,axes[0].unit),
index, integral, self.toXLink() ) )
minProb = min( [xy.toPointwise_withLinearXYs( accuracy = XYsModule.defaultAccuracy, upperEps = 1e-8 ).rangeMin for xy in energy_in] )
if minProb < 0:
warnings.append( warning.negativeProbability( PQU.PQU(energy_in.value,axes[0].unit),
value = minProb, obj = energy_in ) )
if self.interpolationQualifier is standardsModule.interpolation.unitBaseToken:
# check for more than one outgoing distribution integrating to 0 at high end of incident energies
integrals = [eout.integrate() for eout in energy_in]
for idx, integral in enumerate(integrals[::-1]):
if integral != 0: break
if idx > 1:
warnings.append( warning.extraOutgoingEnergy( PQU.PQU(energy_in.value,axes[-1].unit),
obj = energy_in ) )
return warnings
[docs] def normalize( self, insitu = True ) :
n = self
if( not( insitu ) ) : n = self.copy( )
for E_MuEpPs in n :
sum = float( E_MuEpPs.integrate( ) )
for muEpPs in E_MuEpPs : muEpPs.setData( muEpPs / sum )
return( n )
[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']
productFrame = tempInfo['productFrame']
if( verbosity > 2 ) : print '%sGrouping %s' % ( indent, self.moniker )
TM_1, TM_E = transferMatricesModule.EEpMuP_TransferMatrix( style, tempInfo, productFrame, tempInfo['crossSection'], self,
tempInfo['multiplicity'], comment = tempInfo['transferMatrixComment'] + ' outgoing data for %s' % tempInfo['productLabel'] )
return( groupModule.TMs2Form( style, tempInfo, TM_1, TM_E ) )
[docs] def toPointwise_withLinearXYs( self, **kwargs ) :
return( multiD_XYsModule.XYs3d.toPointwise_withLinearXYs( self, cls = XYs3d, **kwargs ) )
[docs] def to_xs_pdf_cdf1d( self, style, tempInfo, indent ) :
from . import energy as energyModule
from . import energyAngularMC as energyAngularMCModule
energy = energyModule.XYs2d( axes = self.axes )
for PofMuGivenEp in self :
data = energyModule.XYs1d( [ [ PofMu.value, float( PofMu.integrate( ) ) ] for PofMu in PofMuGivenEp ], interpolation = PofMuGivenEp.interpolation )
energy.append( energyModule.xs_pdf_cdf1d.fromXYs( data, value = PofMuGivenEp.value ) )
xys3d = energyAngularMCModule.XYs3d( axes = self.axes )
for PofMuGivenEp in self :
xys2d = energyAngularMCModule.XYs2d( value = PofMuGivenEp.value )
for PofMu in PofMuGivenEp :
_PofMu = PofMu
if( isinstance( PofMu, Legendre ) ) :
if( PofMu[0] == 0 ) : _PofMu = XYs1d( [ [ -1, 0.5 ], [ 1, 0.5 ] ] )
_PofMu = _PofMu.toPointwise_withLinearXYs( accuracy = XYsModule.defaultAccuracy, upperEps = 1e-8 )
xys2d.append( energyAngularMCModule.xs_pdf_cdf1d.fromXYs( _PofMu, PofMu.value ) )
xys3d.append( xys2d )
return( energyAngularMCModule.energy( energy ), energyAngularMCModule.energyAngular( xys3d ) )
[docs] @staticmethod
def allowedSubElements( ) :
return( ( XYs2d, ) )
[docs]class regions3d( subform, regionsModule.regions3d ) :
def __init__( self, **kwargs ) :
regionsModule.regions3d.__init__( self, **kwargs )
subform.__init__( self )
[docs] def calculateAverageProductData( self, style, indent = '', **kwargs ) :
raise Exception( 'Not implemented' )
[docs] def check( self, info ) :
raise Exception( 'Not implemented' )
[docs] def normalize( self, insitu = True ) :
n = self
if( not( insitu ) ) : n = self.copy( )
for region in n : region.normalize( insitu = True )
return( n )
[docs] def toPointwise_withLinearXYs( self, **kwargs ) :
raise Exception( 'Not implemented' )
[docs] @staticmethod
def allowedSubElements( ) :
return( ( XYs3d, ) )