Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
MCGIDI_distributions.cc
Go to the documentation of this file.
1/*
2# <<BEGIN-copyright>>
3# Copyright 2019, Lawrence Livermore National Security, LLC.
4# This file is part of the gidiplus package (https://github.com/LLNL/gidiplus).
5# gidiplus is licensed under the MIT license (see https://opensource.org/licenses/MIT).
6# SPDX-License-Identifier: MIT
7# <<END-copyright>>
8*/
9
10#include "MCGIDI.hpp"
11
12namespace MCGIDI {
13
14namespace Distributions {
15
16LUPI_HOST_DEVICE static double coherentPhotoAtomicScatteringIntegrateSub( int a_n, double a_a, double a_logX, double a_energy1, double a_y1, double a_energy2, double a_y2 );
17static LUPI_HOST Distribution *parseGIDI2( GIDI::Distributions::Distribution const &a_GIDI_distribution, SetupInfo &a_setupInfo,
18 Transporting::MC const &a_settings );
19
20/*! \class Distribution
21 * This class is the base class for all distribution forms.
22 */
23
24/* *********************************************************************************************************//**
25 * Default constructor used when broadcasting a Protare as needed by MPI or GPUs.
26 ***********************************************************************************************************/
27
29 m_type( Type::none ),
30 m_productFrame( GIDI::Frame::lab ),
31 m_projectileMass( 0.0 ),
32 m_targetMass( 0.0 ),
33 m_productMass( 0.0 ) {
34
35}
36
37/* *********************************************************************************************************//**
38 * @param a_type [in] The Type of the distribution.
39 * @param a_distribution [in] The GIDI::Distributions::Distribution instance whose data is to be used to construct *this*.
40 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
41 ***********************************************************************************************************/
42
44 m_type( a_type ),
45 m_productFrame( a_distribution.productFrame( ) ),
46 m_projectileMass( a_setupInfo.m_protare.projectileMass( ) ),
47 m_targetMass( a_setupInfo.m_protare.targetMass( ) ),
48 m_productMass( a_setupInfo.m_product1Mass ) { // Includes nuclear excitation energy.
49
50}
51
52/* *********************************************************************************************************//**
53 * @param a_type [in] The Type of the distribution.
54 * @param a_productFrame [in] The frame of the product's data for distribution.
55 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
56 ***********************************************************************************************************/
57
58LUPI_HOST Distribution::Distribution( Type a_type, GIDI::Frame a_productFrame, SetupInfo &a_setupInfo ) :
59 m_type( a_type ),
60 m_productFrame( a_productFrame ),
61 m_projectileMass( a_setupInfo.m_protare.projectileMass( ) ),
62 m_targetMass( a_setupInfo.m_protare.targetMass( ) ),
63 m_productMass( a_setupInfo.m_product1Mass ) { // Includes nuclear excitation energy.
64
65}
66
67/* *********************************************************************************************************//**
68 ***********************************************************************************************************/
69
73
74/* *********************************************************************************************************//**
75 * This method calls the **setModelDBRC_data2* method if the distribution is AngularTwoBody, otherwise it * executes a throw.
76 *
77 * @param a_modelDBRC_data [in] The instance storing data needed to treat the DRRC upscatter mode.
78 ***********************************************************************************************************/
79
81
82 if( type( ) != Type::angularTwoBody ) throw std::runtime_error( "Setting ModelDBRC_data for non-two body distribution is not allowed." );
83
84 static_cast<AngularTwoBody *>( this )->setModelDBRC_data2( a_modelDBRC_data );
85}
86
87/* *********************************************************************************************************//**
88 * This method serializes *this* for broadcasting as needed for MPI and GPUs. The method can count the number of required
89 * bytes, pack *this* or unpack *this* depending on *a_mode*.
90 *
91 * @param a_buffer [in] The buffer to read or write data to depending on *a_mode*.
92 * @param a_mode [in] Specifies the action of this method.
93 ***********************************************************************************************************/
94
96
97 int distributionType = distributionTypeToInt( m_type );
98 DATA_MEMBER_INT( distributionType, a_buffer, a_mode );
99 if( a_mode == LUPI::DataBuffer::Mode::Unpack ) m_type = intToDistributionType( distributionType );
100
101 int frame = 0;
102 if( m_productFrame == GIDI::Frame::centerOfMass ) frame = 1;
103 DATA_MEMBER_INT( frame, a_buffer, a_mode );
104 m_productFrame = GIDI::Frame::lab;
105 if( frame == 1 ) m_productFrame = GIDI::Frame::centerOfMass;
106
107 DATA_MEMBER_DOUBLE( m_projectileMass, a_buffer, a_mode );
108 DATA_MEMBER_DOUBLE( m_targetMass, a_buffer, a_mode );
109 DATA_MEMBER_DOUBLE( m_productMass, a_buffer, a_mode );
110}
111
112/*! \class AngularTwoBody
113 * This class represents the distribution for an outgoing product for a two-body interaction.
114 */
115
116/* *********************************************************************************************************//**
117 * Base contructor.
118 ***********************************************************************************************************/
119
121 m_residualMass( 0.0 ),
122 m_Q( 0.0 ),
123 m_twoBodyThreshold( 0.0 ),
124 m_Upscatter( false ),
125 m_angular( nullptr ),
126 m_modelDBRC_data( nullptr ) {
127
128}
129
130/* *********************************************************************************************************//**
131 * @param a_angularTwoBody [in] The GIDI::Distributions::AngularTwoBody instance whose data is to be used to construct *this*.
132 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
133 ***********************************************************************************************************/
134
136 Distribution( Type::angularTwoBody, a_angularTwoBody, a_setupInfo ),
137 m_residualMass( a_setupInfo.m_product2Mass ), // Includes nuclear excitation energy.
138 m_Q( a_setupInfo.m_Q ),
139 m_twoBodyThreshold( a_setupInfo.m_reaction->twoBodyThreshold( ) ),
140 m_Upscatter( false ),
141 m_angular( Probabilities::parseProbability2d_d1( a_angularTwoBody.angular( ), &a_setupInfo ) ),
142 m_modelDBRC_data( nullptr ) {
143
144 if( a_setupInfo.m_protare.projectileIntid( ) == PoPI::Intids::neutron ) {
145 m_Upscatter = a_setupInfo.m_reaction->ENDF_MT( ) == 2;
146 }
147}
148
149/* *********************************************************************************************************//**
150 ***********************************************************************************************************/
151
153
154 delete m_angular;
155 delete m_modelDBRC_data;
156}
157
158/* *********************************************************************************************************//**
159 * This method serializes *this* for broadcasting as needed for MPI and GPUs. The method can count the number of required
160 * bytes, pack *this* or unpack *this* depending on *a_mode*.
161 *
162 * @param a_buffer [in] The buffer to read or write data to depending on *a_mode*.
163 * @param a_mode [in] Specifies the action of this method.
164 ***********************************************************************************************************/
165
167
168 Distribution::serialize( a_buffer, a_mode );
169 DATA_MEMBER_DOUBLE( m_residualMass, a_buffer, a_mode );
170 DATA_MEMBER_DOUBLE( m_Q, a_buffer, a_mode );
171 DATA_MEMBER_DOUBLE( m_twoBodyThreshold, a_buffer, a_mode );
172 DATA_MEMBER_INT( m_Upscatter, a_buffer, a_mode );
173
174 m_angular = serializeProbability2d_d1( a_buffer, a_mode, m_angular );
175 m_modelDBRC_data = serializeModelDBRC_data( a_buffer, a_mode, m_modelDBRC_data );
176}
177
178/* *********************************************************************************************************//**
179 * This method sets *this* *m_modelDBRC_data* to *a_modelDBRC_data*. It also deletes the current *m_modelDBRC_data* member.
180 *
181 * @param a_modelDBRC_data [in] The instance storing data needed to treat the DRRC upscatter mode.
182 ***********************************************************************************************************/
183
185
186 delete m_modelDBRC_data;
187 m_modelDBRC_data = a_modelDBRC_data;
188}
189
190/*! \class Uncorrelated
191 * This class represents the distribution for an outgoing product for which the distribution is the product of uncorrelated
192 * angular (i.e., P(mu|E)) and energy (i.e., P(E'|E)) distributions.
193 */
194
195/* *********************************************************************************************************//**
196 * Default constructor used when broadcasting a Protare as needed by MPI or GPUs.
197 ***********************************************************************************************************/
198
200 m_angular( nullptr ),
201 m_energy( nullptr ) {
202
203}
204
205/* *********************************************************************************************************//**
206 * @param a_uncorrelated [in] The GIDI::Distributions::Uncorrelated instance whose data is to be used to construct *this*.
207 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
208 ***********************************************************************************************************/
209
211 Distribution( Type::uncorrelated, a_uncorrelated, a_setupInfo ),
212 m_angular( Probabilities::parseProbability2d_d1( a_uncorrelated.angular( ), nullptr ) ),
213 m_energy( Probabilities::parseProbability2d( a_uncorrelated.energy( ), &a_setupInfo ) ) {
214
215}
216
217/* *********************************************************************************************************//**
218 ***********************************************************************************************************/
219
221
222 delete m_angular;
223 delete m_energy;
224}
225
226/* *********************************************************************************************************//**
227 * This method serializes *this* for broadcasting as needed for MPI and GPUs. The method can count the number of required
228 * bytes, pack *this* or unpack *this* depending on *a_mode*.
229 *
230 * @param a_buffer [in] The buffer to read or write data to depending on *a_mode*.
231 * @param a_mode [in] Specifies the action of this method.
232 ***********************************************************************************************************/
233
235
236 Distribution::serialize( a_buffer, a_mode );
237 m_angular = serializeProbability2d_d1( a_buffer, a_mode, m_angular );
238 m_energy = serializeProbability2d( a_buffer, a_mode, m_energy );
239}
240
241/*! \class Branching3d
242 * This class represents the distribution for an outgoing product for which the distribution is the product of uncorrelated
243 * angular (i.e., P(mu|E)) and energy (i.e., P(E'|E)) distributions.
244 */
245
246/* *********************************************************************************************************//**
247 * Default constructor used when broadcasting a Protare as needed by MPI or GPUs.
248 ***********************************************************************************************************/
249
251 m_initialStateIndex( -1 ) {
252
253}
254
255/* *********************************************************************************************************//**
256 * @param a_branching3d [in] The GIDI::Distributions::Branching3dinstance whose data is to be used to construct *this*.
257 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
258 ***********************************************************************************************************/
259
261 Distribution( Type::branching3d, a_branching3d, a_setupInfo ),
262 m_initialStateIndex( -1 ) {
263
264 auto iter = a_setupInfo.m_stateNamesToIndices.find( a_branching3d.initialState( ) );
265 if( iter == a_setupInfo.m_stateNamesToIndices.end( ) ) {
266 std::string message( "Branching3d: initial state not found: pid = '" + a_branching3d.initialState( ) + "'." );
267 throw std::runtime_error( message.c_str( ) );
268 }
269 m_initialStateIndex = iter->second;
270
271 a_setupInfo.m_initialStateIndex = m_initialStateIndex;
272}
273
274/* *********************************************************************************************************//**
275 ***********************************************************************************************************/
276
280
281
282/* *********************************************************************************************************//**
283 * This method serializes *this* for broadcasting as needed for MPI and GPUs. The method can count the number of required
284 * bytes, pack *this* or unpack *this* depending on *a_mode*.
285 *
286 * @param a_buffer [in] The buffer to read or write data to depending on *a_mode*.
287 * @param a_mode [in] Specifies the action of this method.
288 ***********************************************************************************************************/
289
291
292 Distribution::serialize( a_buffer, a_mode );
293 DATA_MEMBER_INT( m_initialStateIndex, a_buffer, a_mode );
294}
295
296/*! \class EnergyAngularMC
297 * This class represents the distribution for an outgoing particle where the distribution is give as
298 * P(E'|E) * P(mu|E,E') where E is the projectile's energy, E' is the product's outgoing energy, mu is the
299 * cosine of the product's outgoing angle relative to the projectile's velocity, P(E'|E) is the probability for E' given E
300 * and (P(mu|E,E') is the probability for mu given E and E'.
301 */
302
303/* *********************************************************************************************************//**
304 * Default constructor used when broadcasting a Protare as needed by MPI or GPUs.
305 ***********************************************************************************************************/
306
308 m_energy( nullptr ),
309 m_angularGivenEnergy( nullptr ) {
310
311}
312
313/* *********************************************************************************************************//**
314 * @param a_energyAngularMC [in] The GIDI::Distributions::EnergyAngularMC instance whose data is to be used to construct *this*.
315 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
316 ***********************************************************************************************************/
317
319 Distribution( Type::energyAngularMC, a_energyAngularMC, a_setupInfo ),
320 m_energy( Probabilities::parseProbability2d_d1( a_energyAngularMC.energy( ), nullptr ) ),
321 m_angularGivenEnergy( Probabilities::parseProbability3d( a_energyAngularMC.energyAngular( ) ) ) {
322
323}
324
325/* *********************************************************************************************************//**
326 ***********************************************************************************************************/
327
329
330 delete m_energy;
331 delete m_angularGivenEnergy;
332}
333
334/* *********************************************************************************************************//**
335 * This method serializes *this* for broadcasting as needed for MPI and GPUs. The method can count the number of required
336 * bytes, pack *this* or unpack *this* depending on *a_mode*.
337 *
338 * @param a_buffer [in] The buffer to read or write data to depending on *a_mode*.
339 * @param a_mode [in] Specifies the action of this method.
340 ***********************************************************************************************************/
341
343
344 Distribution::serialize( a_buffer, a_mode );
345 m_energy = serializeProbability2d_d1( a_buffer, a_mode, m_energy );
346 m_angularGivenEnergy = serializeProbability3d( a_buffer, a_mode, m_angularGivenEnergy );
347}
348
349/*! \class AngularEnergyMC
350 * This class represents the distribution for an outgoing particle where the distribution is give as
351 * P(mu|E) * P(E'|E,mu) where E is the projectile's energy, E' is the product's outgoing energy, mu is the
352 * cosine of the product's outgoing angle relative to the projectile's velocity, P(mu|E) is the probability for mu given E
353 * and (P(E'|E,mu) is the probability for E' given E and mu.
354 */
355
356/* *********************************************************************************************************//**
357 * Default constructor used when broadcasting a Protare as needed by MPI or GPUs.
358 ***********************************************************************************************************/
359
361 m_angular( nullptr ),
362 m_energyGivenAngular( nullptr ) {
363
364}
365
366/* *********************************************************************************************************//**
367 * @param a_angularEnergyMC [in] The GIDI::Distributions::AngularEnergyMC instance whose data is to be used to construct *this*.
368 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
369 ***********************************************************************************************************/
370
372 Distribution( Type::angularEnergyMC, a_angularEnergyMC, a_setupInfo ),
373 m_angular( Probabilities::parseProbability2d_d1( a_angularEnergyMC.angular( ), nullptr ) ),
374 m_energyGivenAngular( Probabilities::parseProbability3d( a_angularEnergyMC.angularEnergy( ) ) ) {
375
376}
377
378/* *********************************************************************************************************//**
379 ***********************************************************************************************************/
380
382
383 delete m_angular;
384 delete m_energyGivenAngular;
385}
386
387/* *********************************************************************************************************//**
388 * This method serializes *this* for broadcasting as needed for MPI and GPUs. The method can count the number of required
389 * bytes, pack *this* or unpack *this* depending on *a_mode*.
390 *
391 * @param a_buffer [in] The buffer to read or write data to depending on *a_mode*.
392 * @param a_mode [in] Specifies the action of this method.
393 ***********************************************************************************************************/
394
396
397 Distribution::serialize( a_buffer, a_mode );
398 m_angular = serializeProbability2d_d1( a_buffer, a_mode, m_angular );
399 m_energyGivenAngular = serializeProbability3d( a_buffer, a_mode, m_energyGivenAngular );
400}
401
402/*! \class KalbachMann
403 * This class represents the distribution for an outgoing product whose distribution is represented by Kalbach-Mann systematics.
404 */
405
406/* *********************************************************************************************************//**
407 * Default constructor used when broadcasting a Protare as needed by MPI or GPUs.
408 ***********************************************************************************************************/
409
411 m_energyToMeVFactor( 0.0 ),
412 m_eb_massFactor( 0.0 ),
413 m_f( nullptr ),
414 m_r( nullptr ),
415 m_a( nullptr ) {
416
417}
418
419/* *********************************************************************************************************//**
420 * @param a_KalbachMann [in] The GIDI::Distributions::KalbachMann instance whose data is to be used to construct *this*.
421 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
422 ***********************************************************************************************************/
423
425 Distribution( Type::KalbachMann, a_KalbachMann, a_setupInfo ),
426 m_energyToMeVFactor( 1 ), // FIXME.
427 m_eb_massFactor( 1 ), // FIXME.
428 m_f( Probabilities::parseProbability2d_d1( a_KalbachMann.f( ), nullptr ) ),
429 m_r( Functions::parseFunction2d( a_KalbachMann.r( ) ) ),
430 m_a( Functions::parseFunction2d( a_KalbachMann.a( ) ) ) {
431
432}
433
434/* *********************************************************************************************************//**
435 ***********************************************************************************************************/
436
438
439 delete m_f;
440 delete m_r;
441 delete m_a;
442}
443
444/* *********************************************************************************************************//**
445 * This method evaluates the Kalbach-Mann formalism at the projectile energy a_energy, and outgoing product energy a_energyOut and a_mu.
446 *
447 * @param a_energy [in] The energy of the projectile in the lab frame.
448 * @param a_energyOut [in] The energy of the product in the center-of-mass frame.
449 * @param a_mu [in] The mu of the product in the center-of-mass frame.
450 ***********************************************************************************************************/
451
452LUPI_HOST_DEVICE double KalbachMann::evaluate( double a_energy, double a_energyOut, double a_mu ) {
453
454// double f_0 = m_f->evaluate( a_energy, a_energyOut );
455 double rValue = m_r->evaluate( a_energy, a_energyOut );
456 double aValue = m_a->evaluate( a_energy, a_energyOut );
457// double pdf_val = aValue * f_0 / 2.0 / sinh( aValue ) * (cosh(aValue * a_mu) + rValue * sinh( aValue * a_mu ) ); // double-differential PDF for a_energyOut and a_mu (Eq. 6.4 in ENDF-102, 2012)
458 double pdf_val = aValue / ( 2.0 * sinh( aValue ) ) * ( cosh( aValue * a_mu ) + rValue * sinh( aValue * a_mu ) ); // double-differential PDF for a_energyOut and mu (Eq. 6.4 in ENDF-102, 2012)
459 return pdf_val;
460}
461
462/* *********************************************************************************************************//**
463 * This method serializes *this* for broadcasting as needed for MPI and GPUs. The method can count the number of required
464 * bytes, pack *this* or unpack *this* depending on *a_mode*.
465 *
466 * @param a_buffer [in] The buffer to read or write data to depending on *a_mode*.
467 * @param a_mode [in] Specifies the action of this method.
468 ***********************************************************************************************************/
469
471
472 Distribution::serialize( a_buffer, a_mode );
473 DATA_MEMBER_DOUBLE( m_energyToMeVFactor, a_buffer, a_mode );
474 DATA_MEMBER_DOUBLE( m_eb_massFactor, a_buffer, a_mode );
475
476 m_f = serializeProbability2d_d1( a_buffer, a_mode, m_f );
477 m_r = serializeFunction2d( a_buffer, a_mode, m_r );
478 m_a = serializeFunction2d( a_buffer, a_mode, m_a );
479}
480
481/*! \class CoherentPhotoAtomicScattering
482 * This class represents the distribution for an outgoing photon via coherent photo-atomic elastic scattering.
483 */
484
485/* *********************************************************************************************************//**
486 * Default constructor used when broadcasting a Protare as needed by MPI or GPUs.
487 ***********************************************************************************************************/
488
490 m_realAnomalousFactor( nullptr ),
491 m_imaginaryAnomalousFactor( nullptr ) {
492
493}
494
495/* *********************************************************************************************************//**
496 * @param a_coherentPhotoAtomicScattering [in] GIDI::Distributions::CoherentPhotoAtomicScattering instance whose data is to be used to construct *this*.
497 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
498 ***********************************************************************************************************/
499
501 Distribution( Type::coherentPhotoAtomicScattering, a_coherentPhotoAtomicScattering, a_setupInfo ),
502 m_anomalousDataPresent( false ),
503 m_realAnomalousFactor( nullptr ),
504 m_imaginaryAnomalousFactor( nullptr ) {
505
506 GUPI::Ancestry const *link = a_coherentPhotoAtomicScattering.findInAncestry( a_coherentPhotoAtomicScattering.href( ) );
509
510 std::string domainUnit;
511 GIDI::Functions::XYs1d const *xys1d0, *xys1d1;
512 std::size_t dataSize = 0, offset = 0;
513
514 GIDI::Functions::Function1dForm const *formFactor = coherentPhotoAtomicScattering.formFactor( );
515 if( formFactor->type( ) == GIDI::FormType::XYs1d ) {
516 xys1d0 = static_cast<GIDI::Functions::XYs1d const *>( formFactor );
517 xys1d1 = xys1d0;
518
519 domainUnit = xys1d0->axes( )[0]->unit( );
520
521 dataSize = xys1d1->size( );
522 offset = 1; }
523 else if( formFactor->type( ) == GIDI::FormType::regions1d ) {
524 GIDI::Functions::Regions1d const *regions1d = static_cast<GIDI::Functions::Regions1d const *>( formFactor );
525 if( regions1d->size( ) != 2 ) throw std::runtime_error( "MCGIDI::CoherentPhotoAtomicScattering::CoherentPhotoAtomicScattering: unsupported form factor size." );
526
527 domainUnit = regions1d->axes( )[0]->unit( );
528
529 GIDI::Functions::Function1dForm const *region0 = (*regions1d)[0];
530 if( region0->type( ) != GIDI::FormType::XYs1d ) throw std::runtime_error( "MCGIDI::CoherentPhotoAtomicScattering::CoherentPhotoAtomicScattering: unsupported form factor for region 0." );
531 xys1d0 = static_cast<GIDI::Functions::XYs1d const *>( region0 );
532 if( xys1d0->size( ) != 2 ) throw std::runtime_error( "MCGIDI::CoherentPhotoAtomicScattering::CoherentPhotoAtomicScattering: unsupported size of region 1 of form factor." );
533
534 GIDI::Functions::Function1dForm const *region1 = (*regions1d)[1];
535 if( region1->type( ) != GIDI::FormType::XYs1d ) throw std::runtime_error( "MCGIDI::CoherentPhotoAtomicScattering::CoherentPhotoAtomicScattering: unsupported form factor for region 1." );
536 xys1d1 = static_cast<GIDI::Functions::XYs1d const *>( region1 );
537
538 dataSize = xys1d1->size( ) + 1; }
539 else {
540 throw std::runtime_error( "MCGIDI::CoherentPhotoAtomicScattering::CoherentPhotoAtomicScattering: unsupported form factor. Must be XYs1d or regions1d." );
541 }
542
543 double domainFactor = 1.0;
544 if( domainUnit == "1/Ang" ) {
545 domainFactor = 0.012398419739640716; } // Converts 'h * c /Ang' to MeV.
546 else if( domainUnit == "1/cm" ) {
547 domainFactor = 0.012398419739640716 * 1e-8; } // Converts 'h * c /cm' to MeV.
548 else {
549 throw std::runtime_error( "MCGIDI::CoherentPhotoAtomicScattering::CoherentPhotoAtomicScattering: unsupported domain unit" );
550 }
551
552 m_energies.resize( dataSize );
553 m_formFactor.resize( dataSize );
554 m_a.resize( dataSize );
555 m_integratedFormFactor.resize( dataSize );
556 m_integratedFormFactorSquared.resize( dataSize );
557 m_probabilityNorm1_1.resize( dataSize );
558 m_probabilityNorm1_3.resize( dataSize );
559 m_probabilityNorm1_5.resize( dataSize );
560 m_probabilityNorm2_1.resize( dataSize );
561 m_probabilityNorm2_3.resize( dataSize );
562 m_probabilityNorm2_5.resize( dataSize );
563
564 std::pair<double, double> xy = (*xys1d0)[0];
565 m_energies[0] = 0.0;
566 m_formFactor[0] = xy.second;
567 m_a[0] = 0.0;
568 m_integratedFormFactor[0] = 0.0;
569 m_integratedFormFactorSquared[0] = 0.0;
570
571 xy = (*xys1d1)[offset];
572 double energy1 = domainFactor * xy.first;
573 double y1 = xy.second;
574 m_energies[1] = energy1;
575 m_formFactor[1] = y1;
576 m_integratedFormFactor[1] = 0.5 * energy1 * energy1 * y1;
577 m_integratedFormFactorSquared[1] = 0.5 * energy1 * energy1 * y1 * y1;
578
579 double sum1 = m_integratedFormFactor[1];
580 double sum2 = m_integratedFormFactorSquared[1];
581 for( std::size_t i1 = 1 + offset; i1 < xys1d1->size( ); ++i1 ) {
582 xy = (*xys1d1)[i1];
583 double energy2 = domainFactor * xy.first;
584 double y2 = xy.second;
585
586 double logEs = log( energy2 / energy1 );
587 double _a = log( y2 / y1 ) / logEs;
588
589 m_energies[i1+1-offset] = energy2;
590 m_formFactor[i1+1-offset] = y2;
591 m_a[i1-offset] = _a;
592
593 sum1 += coherentPhotoAtomicScatteringIntegrateSub( 1, _a, logEs, energy1, y1, energy2, y2 );
594 m_integratedFormFactor[i1+1-offset] = sum1;
595
596 sum2 += coherentPhotoAtomicScatteringIntegrateSub( 1, 2.0 * _a, logEs, energy1, y1 * y1, energy2, y2 * y2 );
597 m_integratedFormFactorSquared[i1+1-offset] = sum2;
598
599 energy1 = energy2;
600 y1 = y2;
601 }
602
603 m_a[m_a.size()-1] = 0.0;
604
605 if( coherentPhotoAtomicScattering.realAnomalousFactor( ) != nullptr ) {
606 m_anomalousDataPresent = true;
607 m_realAnomalousFactor = Functions::parseFunction1d_d1( coherentPhotoAtomicScattering.realAnomalousFactor( ) );
608 m_imaginaryAnomalousFactor = Functions::parseFunction1d_d1( coherentPhotoAtomicScattering.imaginaryAnomalousFactor( ) );
609 }
610
611 m_probabilityNorm1_1[0] = 0.0;
612 m_probabilityNorm1_3[0] = 0.0;
613 m_probabilityNorm1_5[0] = 0.0;
614 m_probabilityNorm2_1[0] = 0.0;
615 m_probabilityNorm2_3[0] = 0.0;
616 m_probabilityNorm2_5[0] = 0.0;
617 energy1 = m_energies[1];
618 y1 = m_formFactor[0];
619 for( std::size_t i1 = 1; i1 < m_probabilityNorm1_1.size( ); ++i1 ) {
620 double energy2 = m_energies[i1];
621 double y2 = m_formFactor[i1];
622 double logEs = log( energy2 / energy1 );
623
624 m_probabilityNorm1_1[i1] = m_probabilityNorm1_1[i1-1] + coherentPhotoAtomicScatteringIntegrateSub( 1, m_a[i1-1], logEs, energy1, y1, energy2, y2 );
625 m_probabilityNorm1_3[i1] = m_probabilityNorm1_3[i1-1] + coherentPhotoAtomicScatteringIntegrateSub( 3, m_a[i1-1], logEs, energy1, y1, energy2, y2 );
626 m_probabilityNorm1_5[i1] = m_probabilityNorm1_5[i1-1] + coherentPhotoAtomicScatteringIntegrateSub( 5, m_a[i1-1], logEs, energy1, y1, energy2, y2 );
627
628 m_probabilityNorm2_1[i1] = m_probabilityNorm2_1[i1-1] + coherentPhotoAtomicScatteringIntegrateSub( 1, 2.0 * m_a[i1-1], logEs, energy1, y1 * y1, energy2, y2 * y2 );
629 m_probabilityNorm2_3[i1] = m_probabilityNorm2_3[i1-1] + coherentPhotoAtomicScatteringIntegrateSub( 3, 2.0 * m_a[i1-1], logEs, energy1, y1 * y1, energy2, y2 * y2 );
630 m_probabilityNorm2_5[i1] = m_probabilityNorm2_5[i1-1] + coherentPhotoAtomicScatteringIntegrateSub( 5, 2.0 * m_a[i1-1], logEs, energy1, y1 * y1, energy2, y2 * y2 );
631
632 energy1 = energy2;
633 y1 = y2;
634 }
635}
636
637/* *********************************************************************************************************//**
638 ***********************************************************************************************************/
639
641
642 delete m_realAnomalousFactor;
643 delete m_imaginaryAnomalousFactor;
644}
645
646/* *********************************************************************************************************//**
647 * This method evaluates the coherent photo-atomic scattering double differentil at the projectile energy a_energy and product cosine of angle a_mu.
648 *
649 * @param a_energyIn [in] The energy of the projectile in the lab frame.
650 * @param a_mu [in] The mu of the product in the center-of-mass frame.
651 ***********************************************************************************************************/
652
653LUPI_HOST_DEVICE double CoherentPhotoAtomicScattering::evaluate( double a_energyIn, double a_mu ) const {
654
655 double probability;
656 int intLowerIndexEnergy = binarySearchVector( a_energyIn, m_energies, true ); // FIXME - need to handle case where lowerIndexEnergy = 0 like in evaluateScatteringFactor.
657 std::size_t lowerIndexEnergy = static_cast<std::size_t>( intLowerIndexEnergy );
658 double _a = m_a[lowerIndexEnergy];
659 double _a_2 = _a * _a;
660 double X1 = m_energies[lowerIndexEnergy];
661 double logEs = log( a_energyIn / X1 );
662 double formFactor_1 = m_formFactor[lowerIndexEnergy];
663 double formFactor_2 = formFactor_1 * formFactor_1;
664 double formFactorEnergyIn_1 = formFactor_1 * pow( a_energyIn / X1, _a );
665 double formFactorEnergyIn_2 = formFactorEnergyIn_1 * formFactorEnergyIn_1;
666 double inverseEnergyIn_1 = 1.0 / a_energyIn;
667 double inverseEnergyIn_2 = inverseEnergyIn_1 * inverseEnergyIn_1;
668 double inverseEnergyIn_3 = inverseEnergyIn_1 * inverseEnergyIn_2;
669 double inverseEnergyIn_4 = inverseEnergyIn_2 * inverseEnergyIn_2;
670 double inverseEnergyIn_5 = inverseEnergyIn_1 * inverseEnergyIn_4;
671 double inverseEnergyIn_6 = inverseEnergyIn_2 * inverseEnergyIn_4;
672
673 double norm = 0.5 * inverseEnergyIn_2 * ( m_probabilityNorm2_1[lowerIndexEnergy] + coherentPhotoAtomicScatteringIntegrateSub( 1, _a_2, logEs, X1, formFactor_2, a_energyIn, formFactorEnergyIn_2 ) )
674 - inverseEnergyIn_4 * ( m_probabilityNorm2_3[lowerIndexEnergy] + coherentPhotoAtomicScatteringIntegrateSub( 3, _a_2, logEs, X1, formFactor_2, a_energyIn, formFactorEnergyIn_2 ) )
675 + inverseEnergyIn_6 * ( m_probabilityNorm2_5[lowerIndexEnergy] + coherentPhotoAtomicScatteringIntegrateSub( 5, _a_2, logEs, X1, formFactor_2, a_energyIn, formFactorEnergyIn_2 ) );
676
677 double realAnomalousFactor = 0.0;
678 double imaginaryAnomalousFactor = 0.0;
679 if( m_anomalousDataPresent ) {
680 realAnomalousFactor = m_realAnomalousFactor->evaluate( a_energyIn );
681 imaginaryAnomalousFactor = m_imaginaryAnomalousFactor->evaluate( a_energyIn );
682 norm += realAnomalousFactor * ( inverseEnergyIn_1 * ( m_probabilityNorm1_1[lowerIndexEnergy] + coherentPhotoAtomicScatteringIntegrateSub( 1, _a, logEs, X1, formFactor_1, a_energyIn, formFactorEnergyIn_1 ) )
683 - 2.0 * inverseEnergyIn_3 * ( m_probabilityNorm1_3[lowerIndexEnergy] + coherentPhotoAtomicScatteringIntegrateSub( 3, _a, logEs, X1, formFactor_1, a_energyIn, formFactorEnergyIn_1 ) )
684 + 2.0 * inverseEnergyIn_5 * ( m_probabilityNorm1_5[lowerIndexEnergy] + coherentPhotoAtomicScatteringIntegrateSub( 5, _a, logEs, X1, formFactor_1, a_energyIn, formFactorEnergyIn_1 ) ) );
685 }
686 norm *= 16.0;
687 norm += 8.0 / 3.0 * ( realAnomalousFactor * realAnomalousFactor + imaginaryAnomalousFactor * imaginaryAnomalousFactor );
688
689 double _formFactor = evaluateFormFactor( a_energyIn, a_mu );
690 probability = ( 1.0 + a_mu * a_mu ) * ( ( _formFactor + realAnomalousFactor ) * ( _formFactor + realAnomalousFactor ) + imaginaryAnomalousFactor * imaginaryAnomalousFactor ) / norm;
691
692 return( probability );
693}
694
695/* *********************************************************************************************************//**
696 * This method evaluates the coherent photo-atomic form factor at the projectile energy a_energy and product cosine of angle a_mu.
697 *
698 * @param a_energyIn [in] The energy of the projectile in the lab frame.
699 * @param a_mu [in] The mu of the product in the center-of-mass frame.
700 ***********************************************************************************************************/
701
702LUPI_HOST_DEVICE double CoherentPhotoAtomicScattering::evaluateFormFactor( double a_energyIn, double a_mu ) const {
703
704 double X = a_energyIn * sqrt( 0.5 * ( 1 - a_mu ) );
705 int intLowerIndex = binarySearchVector( X, m_energies );
706
707 if( intLowerIndex < 1 ) {
708 if( intLowerIndex == 0 ) return( m_formFactor[0] );
709 if( intLowerIndex == -2 ) return( m_formFactor[0] ); // This should never happend for proper a_energyIn and a_mu.
710 return( m_formFactor.back( ) );
711 }
712
713 std::size_t lowerIndex = static_cast<std::size_t>( intLowerIndex );
714
715 return( m_formFactor[lowerIndex] * pow( X / m_energies[lowerIndex] , m_a[lowerIndex] ) );
716}
717
718/* *********************************************************************************************************//**
719 * FIX ME.
720 *
721 * @param a_Z [in]
722 * @param a_a [in]
723 ***********************************************************************************************************/
724
725LUPI_HOST_DEVICE double CoherentPhotoAtomicScattering::Z_a( double a_Z, double a_a ) const {
726
727 if( fabs( a_a ) < 1e-3 ) {
728 double logZ = log( a_Z );
729 double a_logZ = a_a * logZ;
730 return( logZ * ( 1.0 + 0.5 * a_logZ * ( 1.0 + a_logZ / 3.0 * ( 1.0 + 0.25 * a_logZ ) ) ) );
731 }
732 return( ( pow( a_Z, a_a ) - 1.0 ) / a_a );
733}
734
735/* *********************************************************************************************************//**
736 * FIX ME.
737 *
738 * @param a_n [in]
739 * @param a_a [in]
740 * @param a_logX [in]
741 * @param a_energy1 [in]
742 * @param a_y1 [in]
743 * @param a_energy2 [in]
744 * @param a_y2 [in]
745 ***********************************************************************************************************/
746
747LUPI_HOST_DEVICE static double coherentPhotoAtomicScatteringIntegrateSub( int a_n, double a_a, double a_logX, double a_energy1, double a_y1, double a_energy2, double a_y2 ) {
748
749 double epsilon = a_a + a_n + 1.0;
750 double integral = 0.0;
751
752 if( fabs( epsilon ) < 1e-3 ) {
753 double epsilon_logX = epsilon * a_logX;
754 integral = a_y1 * pow( a_energy1, a_n + 1.0 ) * a_logX * ( 1.0 + 0.5 * epsilon_logX * ( 1.0 + epsilon_logX / 3.0 * ( 1.0 + 0.25 * epsilon_logX ) ) ); }
755 else {
756 integral = ( a_y2 * pow( a_energy2, a_n + 1.0 ) - a_y1 * pow( a_energy1, a_n + 1.0 ) ) / epsilon;
757 }
758
759 return( integral );
760}
761
762/* *********************************************************************************************************//**
763 * This method serializes *this* for broadcasting as needed for MPI and GPUs. The method can count the number of required
764 * bytes, pack *this* or unpack *this* depending on *a_mode*.
765 *
766 * @param a_buffer [in] The buffer to read or write data to depending on *a_mode*.
767 * @param a_mode [in] Specifies the action of this method.
768 ***********************************************************************************************************/
769
771
772 Distribution::serialize( a_buffer, a_mode );
773
774 DATA_MEMBER_INT( m_anomalousDataPresent, a_buffer, a_mode );
775 DATA_MEMBER_VECTOR_DOUBLE( m_energies, a_buffer, a_mode );
776 DATA_MEMBER_VECTOR_DOUBLE( m_formFactor, a_buffer, a_mode );
777 DATA_MEMBER_VECTOR_DOUBLE( m_a, a_buffer, a_mode );
778 DATA_MEMBER_VECTOR_DOUBLE( m_integratedFormFactor, a_buffer, a_mode );
779 DATA_MEMBER_VECTOR_DOUBLE( m_integratedFormFactorSquared, a_buffer, a_mode );
780 DATA_MEMBER_VECTOR_DOUBLE( m_probabilityNorm1_1, a_buffer, a_mode );
781 DATA_MEMBER_VECTOR_DOUBLE( m_probabilityNorm1_3, a_buffer, a_mode );
782 DATA_MEMBER_VECTOR_DOUBLE( m_probabilityNorm1_5, a_buffer, a_mode );
783 DATA_MEMBER_VECTOR_DOUBLE( m_probabilityNorm2_1, a_buffer, a_mode );
784 DATA_MEMBER_VECTOR_DOUBLE( m_probabilityNorm2_3, a_buffer, a_mode );
785 DATA_MEMBER_VECTOR_DOUBLE( m_probabilityNorm2_5, a_buffer, a_mode );
786
787 if( m_anomalousDataPresent ) {
788 m_realAnomalousFactor = serializeFunction1d_d1( a_buffer, a_mode, m_realAnomalousFactor );
789 m_imaginaryAnomalousFactor = serializeFunction1d_d1( a_buffer, a_mode, m_imaginaryAnomalousFactor );
790 }
791}
792
793/*! \class IncoherentPhotoAtomicScattering
794 * This class represents the distribution for an outgoing photon via incoherent photo-atomic elastic scattering.
795 */
796
797/* *********************************************************************************************************//**
798 * Default constructor used when broadcasting a Protare as needed by MPI or GPUs.
799 ***********************************************************************************************************/
800
804
805/* *********************************************************************************************************//**
806 * @param a_incoherentPhotoAtomicScattering [in] The GIDI::Distributions::IncoherentPhotoAtomicScattering instance whose data is to be used to construct *this*.
807 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
808 ***********************************************************************************************************/
809
811 SetupInfo &a_setupInfo ) :
812 Distribution( Type::incoherentPhotoAtomicScattering, a_incoherentPhotoAtomicScattering, a_setupInfo ) {
813
814 GUPI::Ancestry const *link = a_incoherentPhotoAtomicScattering.findInAncestry( a_incoherentPhotoAtomicScattering.href( ) );
817
818 std::string domainUnit;
819 GIDI::Functions::XYs1d const *xys1d0, *xys1d1;
820 std::size_t dataSize = 0, offset = 0;
821
822 GIDI::Functions::Function1dForm const *scatteringFactor = incoherentPhotoAtomicScattering.scatteringFactor( );
823 if( scatteringFactor->type( ) == GIDI::FormType::XYs1d ) {
824 xys1d0 = static_cast<GIDI::Functions::XYs1d const *>( scatteringFactor );
825 xys1d1 = xys1d0;
826
827 domainUnit = xys1d0->axes( )[0]->unit( );
828
829 dataSize = xys1d1->size( );
830 offset = 1; }
831 else if( scatteringFactor->type( ) == GIDI::FormType::regions1d ) {
832 GIDI::Functions::Regions1d const *regions1d = static_cast<GIDI::Functions::Regions1d const *>( scatteringFactor );
833 if( regions1d->size( ) != 2 ) throw std::runtime_error( "MCGIDI::CoherentPhotoAtomicScattering::CoherentPhotoAtomicScattering: unsupported form factor size." );
834
835 domainUnit = regions1d->axes( )[0]->unit( );
836
837 GIDI::Functions::Function1dForm const *region0 = (*regions1d)[0];
838 if( region0->type( ) != GIDI::FormType::XYs1d ) throw std::runtime_error( "MCGIDI::CoherentPhotoAtomicScattering::CoherentPhotoAtomicScattering: unsupported form factor for region 0." );
839 xys1d0 = static_cast<GIDI::Functions::XYs1d const *>( region0 );
840 if( xys1d0->size( ) != 2 ) throw std::runtime_error( "MCGIDI::CoherentPhotoAtomicScattering::CoherentPhotoAtomicScattering: unsupported size of region 1 of form factor." );
841
842 GIDI::Functions::Function1dForm const *region1 = (*regions1d)[1];
843 if( region1->type( ) != GIDI::FormType::XYs1d ) throw std::runtime_error( "MCGIDI::CoherentPhotoAtomicScattering::CoherentPhotoAtomicScattering: unsupported form factor for region 1." );
844 xys1d1 = static_cast<GIDI::Functions::XYs1d const *>( region1 );
845
846 dataSize = xys1d1->size( ) + 1; }
847 else {
848 throw std::runtime_error( "MCGIDI::CoherentPhotoAtomicScattering::CoherentPhotoAtomicScattering: unsupported form factor. Must be XYs1d or regions1d." );
849 }
850
851 double domainFactor = 1.0;
852 if( domainUnit == "1/Ang" ) {
853 domainFactor = 0.012398419739640716; } // Converts 'h * c /Ang' to MeV.
854 else if( domainUnit == "1/cm" ) {
855 domainFactor = 0.012398419739640716 * 1e-8; } // Converts 'h * c /cm' to MeV.
856 else {
857 throw std::runtime_error( "MCGIDI::IncoherentPhotoAtomicScattering::IncoherentPhotoAtomicScattering: unsupported domain unit" );
858 }
859
860 m_energies.resize( dataSize );
861 m_scatteringFactor.resize( dataSize );
862 m_a.resize( dataSize );
863
864 std::pair<double, double> xy = (*xys1d0)[0];
865 m_energies[0] = domainFactor * xy.first;
866 m_scatteringFactor[0] = xy.second;
867 m_a[0] = 1.0;
868
869 xy = (*xys1d1)[offset];
870 double energy1 = domainFactor * xy.first;
871 double y1 = xy.second;
872
873 m_energies[1] = energy1;
874 m_scatteringFactor[1] = y1;
875
876 for( std::size_t i1 = 1 + offset; i1 < xys1d1->size( ); ++i1 ) {
877 xy = (*xys1d1)[i1];
878 double energy2 = domainFactor * xy.first;
879 double y2 = xy.second;
880
881 m_energies[i1+1-offset] = energy2;
882 m_scatteringFactor[i1+1-offset] = y2;
883
884 double _a = log( y2 / y1 ) / log( energy2 / energy1 );
885 m_a[i1-offset] = _a;
886
887 energy1 = energy2;
888 y1 = y2;
889 }
890 m_a[m_a.size()-1] = 0.0;
891}
892
893/* *********************************************************************************************************//**
894 ***********************************************************************************************************/
895
899
900/* *********************************************************************************************************//**
901 * FIX ME.
902 *
903 * @param a_energyIn [in]
904 * @param a_mu [in]
905 ***********************************************************************************************************/
906
907LUPI_HOST_DEVICE double IncoherentPhotoAtomicScattering::energyRatio( double a_energyIn, double a_mu ) const {
908
909 double relativeEnergy = a_energyIn / PoPI_electronMass_MeV_c2;
910
911 return( 1.0 / ( 1.0 + relativeEnergy * ( 1.0 - a_mu ) ) );
912}
913
914/* *********************************************************************************************************//**
915 * This method evaluates the Klein-Nishina. FIX ME. This should be a function as it does not use member data.
916 *
917 * @param a_energyIn [in]
918 * @param a_mu [in]
919 ***********************************************************************************************************/
920
921LUPI_HOST_DEVICE double IncoherentPhotoAtomicScattering::evaluateKleinNishina( double a_energyIn, double a_mu ) const {
922
923 double relativeEnergy = a_energyIn / PoPI_electronMass_MeV_c2;
924 double _energyRatio = energyRatio( a_energyIn, a_mu );
925 double one_minus_mu = 1.0 - a_mu;
926
927 double norm = ( 1.0 + 2.0 * relativeEnergy );
928 norm = 2.0 * relativeEnergy * ( 2.0 + relativeEnergy * ( 1.0 + relativeEnergy ) * ( 8.0 + relativeEnergy ) ) / ( norm * norm );
929 norm += ( ( relativeEnergy - 2.0 ) * relativeEnergy - 2.0 ) * log( 1.0 + 2.0 * relativeEnergy );
930 norm /= relativeEnergy * relativeEnergy * relativeEnergy;
931
932 return( _energyRatio * _energyRatio * ( _energyRatio + a_mu * a_mu + relativeEnergy * one_minus_mu * one_minus_mu ) / norm );
933}
934
935/* *********************************************************************************************************//**
936 * This method evaluates the Klein-Nishina.
937 *
938 * @param a_energyIn [in] The energy of the projectile.
939 ***********************************************************************************************************/
940
942
943 int intLowerIndex = binarySearchVector( a_energyIn, m_energies );
944
945 if( intLowerIndex < 1 ) {
946 if( intLowerIndex == -1 ) return( m_scatteringFactor.back( ) );
947 return( m_scatteringFactor[1] * a_energyIn / m_energies[1] );
948 }
949
950 std::size_t lowerIndex = static_cast<std::size_t>( intLowerIndex );
951 return( m_scatteringFactor[lowerIndex] * pow( a_energyIn / m_energies[lowerIndex], m_a[lowerIndex] ) );
952}
953
954/* *********************************************************************************************************//**
955 * This method serializes *this* for broadcasting as needed for MPI and GPUs. The method can count the number of required
956 * bytes, pack *this* or unpack *this* depending on *a_mode*.
957 *
958 * @param a_buffer [in] The buffer to read or write data to depending on *a_mode*.
959 * @param a_mode [in] Specifies the action of this method.
960 ***********************************************************************************************************/
961
963
964 Distribution::serialize( a_buffer, a_mode );
965
966 DATA_MEMBER_VECTOR_DOUBLE( m_energies, a_buffer, a_mode );
967 DATA_MEMBER_VECTOR_DOUBLE( m_scatteringFactor, a_buffer, a_mode );
968 DATA_MEMBER_VECTOR_DOUBLE( m_a, a_buffer, a_mode );
969}
970
971/*! \class IncoherentBoundToFreePhotoAtomicScattering
972 * This class represents the distribution for an outgoing photon via incoherent photo-atomic elastic scattering.
973 */
974
975/* *********************************************************************************************************//**
976 * Default constructor used when broadcasting a Protare as needed by MPI or GPUs.
977 ***********************************************************************************************************/
978
982
983/* *********************************************************************************************************//**
984 * @param a_incoherentBoundToFreePhotoAtomicScattering [in] The GIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering instance whose data is to be used to construct *this*.
985 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
986 ***********************************************************************************************************/
987
989 GIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering const &a_incoherentBoundToFreePhotoAtomicScattering,
990 SetupInfo &a_setupInfo ) :
991 Distribution( Type::incoherentBoundToFreePhotoAtomicScattering, a_incoherentBoundToFreePhotoAtomicScattering, a_setupInfo ),
992 m_bindingEnergy( 0.0 ) {
993
994 GIDI::ProtareSingle const &GIDI_protare = a_setupInfo.m_GIDI_protare;
995 auto monikers = GIDI_protare.styles( ).findAllOfMoniker( GIDI_MonteCarlo_cdfStyleChars );
996 std::string MonteCarlo_cdf = "";
997 if( monikers.size( ) == 1 ) {
998 MonteCarlo_cdf = monikers[0][0]->label( );
999 }
1000
1001 std::string Compton_href = a_incoherentBoundToFreePhotoAtomicScattering.href( );
1002
1003 if( Compton_href.find( MonteCarlo_cdf ) != std::string::npos ) {
1004 const GUPI::Ancestry *link = a_incoherentBoundToFreePhotoAtomicScattering.findInAncestry( Compton_href );
1005 GIDI::DoubleDifferentialCrossSection::IncoherentBoundToFreePhotoAtomicScattering const &dd = *static_cast<GIDI::DoubleDifferentialCrossSection::IncoherentBoundToFreePhotoAtomicScattering const *>( link );
1006 GIDI::Functions::Xs_pdf_cdf1d const *xpcCompton;
1007 GIDI::Functions::Function1dForm const *ComptonProfile = dd.ComptonProfile( );
1008 xpcCompton = static_cast<GIDI::Functions::Xs_pdf_cdf1d const *>( ComptonProfile );
1009
1010 std::vector<double> occupationNumbers = xpcCompton->cdf( );
1011 std::vector<double> pz_grid = xpcCompton->Xs( );
1012 std::size_t dataSize = pz_grid.size( );
1013 m_occupationNumber.resize( dataSize );
1014 m_pz.resize( dataSize );
1015 for( std::size_t index = 0; index < occupationNumbers.size( ); ++index ) {
1016 m_occupationNumber[index] = occupationNumbers[index];
1017 m_pz[index] = pz_grid[index];
1018 }
1019
1020 m_bindingEnergy = a_setupInfo.m_Q;
1021 }
1022}
1023
1024/* *********************************************************************************************************//**
1025 ***********************************************************************************************************/
1026
1030
1031/* *********************************************************************************************************//**
1032 * FIX ME.
1033 *
1034 * @param a_energyIn [in]
1035 * @param a_mu [in]
1036 ***********************************************************************************************************/
1037
1039
1040 double relativeEnergy = a_energyIn / PoPI_electronMass_MeV_c2;
1041
1042 return( 1.0 / ( 1.0 + relativeEnergy * ( 1.0 - a_mu ) ) );
1043}
1044
1045/* *********************************************************************************************************//**
1046 * This method evaluates the Klein-Nishina. FIX ME. This should be a function as it does not use member data.
1047 *
1048 * @param a_energyIn [in]
1049 * @param a_mu [in]
1050 ***********************************************************************************************************/
1051
1053
1054 double relativeEnergy = a_energyIn / PoPI_electronMass_MeV_c2;
1055 double _energyRatio = energyRatio( a_energyIn, a_mu );
1056 double one_minus_mu = 1.0 - a_mu;
1057
1058 double norm = ( 1.0 + 2.0 * relativeEnergy );
1059 norm = 2.0 * relativeEnergy * ( 2.0 + relativeEnergy * ( 1.0 + relativeEnergy ) * ( 8.0 + relativeEnergy ) ) / ( norm * norm );
1060 norm += ( ( relativeEnergy - 2.0 ) * relativeEnergy - 2.0 ) * log( 1.0 + 2.0 * relativeEnergy );
1061 norm /= relativeEnergy * relativeEnergy * relativeEnergy;
1062
1063 return( _energyRatio * _energyRatio * ( _energyRatio + a_mu * a_mu + relativeEnergy * one_minus_mu * one_minus_mu ) / norm );
1064}
1065
1066
1067/* *********************************************************************************************************//**
1068 * This method evaluates the occupation number cdf.
1069 *
1070 * @param a_energyIn [in]
1071 * @param a_mu [in]
1072 ***********************************************************************************************************/
1073
1075
1076 const double alpha_in = a_energyIn / PoPI_electronMass_MeV_c2;
1077 //const double mec = 2.7309245307378233e-22; // m_e * c in SI units
1078 const double alpha_binding = -m_bindingEnergy/PoPI_electronMass_MeV_c2; // BE [MeV] / 0.511 [MeV]
1079 const double pzmax = ( -alpha_binding + alpha_in*(alpha_in - alpha_binding)*(1-a_mu) )/( sqrt( 2*alpha_in*(alpha_in-alpha_binding)*(1-a_mu) + alpha_binding*alpha_binding ) ); // *mec
1080
1081 int intLowerIndex = binarySearchVector( pzmax, m_pz );
1082 std::size_t lowerIndex = static_cast<std::size_t>( intLowerIndex );
1083 int size1 = static_cast<int>( m_occupationNumber.size( ) );
1084
1085 if( intLowerIndex == -1 || intLowerIndex == ( size1 - 1 ) ) {
1086 return( m_occupationNumber.back( ) );
1087 }
1088 if( intLowerIndex == -2 ){
1089 return( m_occupationNumber[0] );
1090 }
1091
1092 return(m_occupationNumber[lowerIndex] + (pzmax-m_pz[lowerIndex])*(m_occupationNumber[lowerIndex+1]-m_occupationNumber[lowerIndex])/(m_pz[lowerIndex+1]-m_pz[lowerIndex]) );
1093
1094}
1095
1096/* *********************************************************************************************************//**
1097 * This method serializes *this* for broadcasting as needed for MPI and GPUs. The method can count the number of required
1098 * bytes, pack *this* or unpack *this* depending on *a_mode*.
1099 *
1100 * @param a_buffer [in] The buffer to read or write data to depending on *a_mode*.
1101 * @param a_mode [in] Specifies the action of this method.
1102 ***********************************************************************************************************/
1103
1105
1106 Distribution::serialize( a_buffer, a_mode );
1107
1108 DATA_MEMBER_VECTOR_DOUBLE( m_pz, a_buffer, a_mode );
1109 DATA_MEMBER_VECTOR_DOUBLE( m_occupationNumber, a_buffer, a_mode );
1110
1111}
1112
1113/*
1114======================================================================================================
1115========== IncoherentPhotoAtomicScatteringElectron ==========
1116======================================================================================================
1117*/
1118
1119/*! \class IncoherentPhotoAtomicScatteringElectron
1120 * This class represents the distribution for the outgoing electron for incoherent photo-atomic scattering.
1121 */
1122
1123/* *********************************************************************************************************//**
1124 * Plain constructor.
1125 ***********************************************************************************************************/
1126
1130
1131/* *********************************************************************************************************//**
1132 * Constructor.
1133 *
1134 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
1135 ***********************************************************************************************************/
1136
1141
1142/* *********************************************************************************************************//**
1143 * Destructor.
1144 ***********************************************************************************************************/
1145
1149
1150/* *********************************************************************************************************//**
1151 * This method serializes *this* for broadcasting as needed for MPI and GPUs. The method can count the number of required
1152 * bytes, pack *this* or unpack *this* depending on *a_mode*.
1153 *
1154 * @param a_buffer [in] The buffer to read or write data to depending on *a_mode*.
1155 * @param a_mode [in] Specifies the action of this method.
1156 ***********************************************************************************************************/
1157
1162
1163/*
1164======================================================================================================
1165========== PairProductionGamma ==========
1166======================================================================================================
1167*/
1168
1169/*! \class PairProductionGamma
1170 * This class represents the distribution for an outgoing photon the is the result of an electron annihilating with a positron.
1171 */
1172
1173/* *********************************************************************************************************//**
1174 * Basic constructor.
1175 ***********************************************************************************************************/
1176
1178 m_firstSampled( false ) {
1179
1180}
1181
1182/* *********************************************************************************************************//**
1183 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
1184 * @param a_firstSampled [in] FIX ME
1185 ***********************************************************************************************************/
1186
1188 Distribution( Type::pairProductionGamma, GIDI::Frame::lab, a_setupInfo ),
1189 m_firstSampled( a_firstSampled ) {
1190
1191}
1192
1193/* *********************************************************************************************************//**
1194 ***********************************************************************************************************/
1195
1199
1200
1201/* *********************************************************************************************************//**
1202 * This method serializes *this* for broadcasting as needed for MPI and GPUs. The method can count the number of required
1203 * bytes, pack *this* or unpack *this* depending on *a_mode*.
1204 *
1205 * @param a_buffer [in] The buffer to read or write data to depending on *a_mode*.
1206 * @param a_mode [in] Specifies the action of this method.
1207 ***********************************************************************************************************/
1208
1210
1211 Distribution::serialize( a_buffer, a_mode );
1212
1213 DATA_MEMBER_INT( m_firstSampled, a_buffer, a_mode );
1214}
1215
1216/*! \class CoherentElasticTNSL
1217 * This class represents the distribution for an outgoing product whose distribution is TNSL coherent elastic scattering.
1218 * This class samples directly from the Debye/Waller function.
1219 */
1220
1221/* *********************************************************************************************************//**
1222 * Constructor for the CoherentElasticTNSL class.
1223 ***********************************************************************************************************/
1224
1226 m_temperatureInterpolation( Interpolation::LINLIN ) {
1227
1228}
1229
1230/* *********************************************************************************************************//**
1231 * Constructor for the CoherentElasticTNSL class.
1232 *
1233 * @param a_coherentElasticTNSL [in] GIDI::CoherentElastic instance containing the Debye/Waller data.
1234 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
1235 ***********************************************************************************************************/
1236
1238 SetupInfo &a_setupInfo ) :
1239 Distribution( Type::coherentElasticTNSL, GIDI::Frame::lab, a_setupInfo ),
1240 m_temperatureInterpolation( Interpolation::LINLIN ) {
1241
1243 GIDI::Functions::Gridded2d const *gridded2d = dynamic_cast<GIDI::Functions::Gridded2d const *>( s_table.function2d( ) );
1244 GIDI::Axes const &axes = gridded2d->axes( );
1245
1246 GIDI::Grid const *axis = dynamic_cast<GIDI::Grid const *>( axes[0] );
1247 m_temperatureInterpolation = GIDI2MCGIDI_interpolation( ptwXY_stringToInterpolation( axis->interpolation( ).c_str( ) ) );
1248
1249 double temperatureToMeV_K = 1.0;
1250 if( axis->unit( ) == "K" ) temperatureToMeV_K = 8.617330337217212e-11; // This is a kludge until units are properly supported.
1251 nf_Buffer<double> const *grid = &axis->values( );
1252 m_temperatures.resize( grid->size( ) );
1253 for( std::size_t index = 0; index < grid->size( ); ++index ) m_temperatures[index] = temperatureToMeV_K * (*grid)[index];
1254
1255 axis = dynamic_cast<GIDI::Grid const *>( axes[1] );
1256 grid = &axis->values( );
1257 m_energies.resize( grid->size( ) );
1258 for( std::size_t index = 0; index < grid->size( ); ++index ) m_energies[index] = (*grid)[index];
1259
1260 GIDI::Array::FullArray fullArray = gridded2d->array( ).constructArray( );
1261 m_S_table.resize( fullArray.size( ) );
1262 for( std::size_t index = 0; index < fullArray.m_flattenedValues.size( ); ++index ) m_S_table[index] = fullArray.m_flattenedValues[index];
1263}
1264
1265/* *********************************************************************************************************//**
1266 * This method serializes *this* for broadcasting as needed for MPI and GPUs. The method can count the number of required
1267 * bytes, pack *this* or unpack *this* depending on *a_mode*.
1268 *
1269 * @param a_buffer [in] The buffer to read or write data to depending on *a_mode*.
1270 * @param a_mode [in] Specifies the action of this method.
1271 ***********************************************************************************************************/
1272
1274
1275 Distribution::serialize( a_buffer, a_mode );
1276
1277 DATA_MEMBER_VECTOR_DOUBLE( m_temperatures, a_buffer, a_mode );
1278 DATA_MEMBER_VECTOR_DOUBLE( m_energies, a_buffer, a_mode );
1279 DATA_MEMBER_VECTOR_DOUBLE( m_S_table, a_buffer, a_mode );
1280
1281 int interpolation = 0;
1282 if( a_mode != LUPI::DataBuffer::Mode::Unpack ) {
1283 switch( m_temperatureInterpolation ) {
1284 case Interpolation::FLAT :
1285 break;
1287 interpolation = 1;
1288 break;
1290 interpolation = 2;
1291 break;
1293 interpolation = 3;
1294 break;
1296 interpolation = 4;
1297 break;
1299 interpolation = 5;
1300 break;
1301 }
1302 }
1303 DATA_MEMBER_INT( interpolation, a_buffer, a_mode );
1304 if( a_mode == LUPI::DataBuffer::Mode::Unpack ) {
1305 switch( interpolation ) {
1306 case 0 :
1307 m_temperatureInterpolation = Interpolation::FLAT;
1308 break;
1309 case 1 :
1310 m_temperatureInterpolation = Interpolation::LINLIN;
1311 break;
1312 case 2 :
1313 m_temperatureInterpolation = Interpolation::LINLOG;
1314 break;
1315 case 3 :
1316 m_temperatureInterpolation = Interpolation::LOGLIN;
1317 break;
1318 case 4 :
1319 m_temperatureInterpolation = Interpolation::LOGLOG;
1320 break;
1321 case 5 :
1322 m_temperatureInterpolation = Interpolation::OTHER;
1323 break;
1324 }
1325 }
1326}
1327
1328/*! \class IncoherentElasticTNSL
1329 * This class represents the distribution for an outgoing product whose distribution is TNSL incoherent elastic scattering.
1330 * This class samples directly from the Debye/Waller function.
1331 */
1332
1333/* *********************************************************************************************************//**
1334 * Constructor for the IncoherentElasticTNSL class.
1335 ***********************************************************************************************************/
1336
1338 m_temperatureToMeV_K( 1.0 ),
1339 m_DebyeWallerIntegral( nullptr ) {
1340
1341}
1342
1343/* *********************************************************************************************************//**
1344 * Constructor for the IncoherentElasticTNSL class.
1345 *
1346 * @param a_incoherentElasticTNSL [in] GIDI::IncoherentElastic instance containing the Debye/Waller data.
1347 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
1348 ***********************************************************************************************************/
1349
1351 SetupInfo &a_setupInfo ) :
1352 Distribution( Type::incoherentElasticTNSL, GIDI::Frame::lab, a_setupInfo ),
1353 m_temperatureToMeV_K( 1.0 ),
1354 m_DebyeWallerIntegral( nullptr ) {
1355
1357 m_DebyeWallerIntegral = Functions::parseFunction1d_d1( debyeWallerIntegral.function1d( ) );
1358 GIDI::Axes const &axes = debyeWallerIntegral.function1d( )->axes( );
1359 GIDI::Axis const *axis = dynamic_cast<GIDI::Axis const *>( axes[0] );
1360 if( axis->unit( ) == "K" ) m_temperatureToMeV_K = 8.617330337217212e-11; // This is a kludge until units are properly supported.
1361}
1362
1363/* *********************************************************************************************************//**
1364 * This method serializes *this* for broadcasting as needed for MPI and GPUs. The method can count the number of required
1365 * bytes, pack *this* or unpack *this* depending on *a_mode*.
1366 *
1367 * @param a_buffer [in] The buffer to read or write data to depending on *a_mode*.
1368 * @param a_mode [in] Specifies the action of this method.
1369 ***********************************************************************************************************/
1370
1372
1373 Distribution::serialize( a_buffer, a_mode );
1374
1375 DATA_MEMBER_DOUBLE( m_temperatureToMeV_K, a_buffer, a_mode );
1376 m_DebyeWallerIntegral = serializeFunction1d_d1( a_buffer, a_mode, m_DebyeWallerIntegral );
1377}
1378
1379/*! \class Unspecified
1380 * This class represents the distribution for an outgoing product whose distribution is not specified.
1381 */
1382
1386
1387/* *********************************************************************************************************//**
1388 * @param a_distribution [in] The GIDI::Distributions::Unspecified instance whose data is to be used to construct *this*.
1389 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
1390 ***********************************************************************************************************/
1391
1393 Distribution( Type::unspecified, a_distribution, a_setupInfo ) {
1394
1395}
1396
1397/* *********************************************************************************************************//**
1398 ***********************************************************************************************************/
1399
1403
1404/* *********************************************************************************************************//**
1405 * This method serializes *this* for broadcasting as needed for MPI and GPUs. The method can count the number of required
1406 * bytes, pack *this* or unpack *this* depending on *a_mode*.
1407 *
1408 * @param a_buffer [in] The buffer to read or write data to depending on *a_mode*.
1409 * @param a_mode [in] Specifies the action of this method.
1410 ***********************************************************************************************************/
1411
1413
1414 Distribution::serialize( a_buffer, a_mode );
1415}
1416
1417
1418/* *********************************************************************************************************//**
1419 * This function is used to call the proper distribution constructor for *a_distribution*.
1420 *
1421 * @param a_distribution [in] The GIDI::Protare whose data is to be used to construct *this*.
1422 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
1423 * @param a_settings [in] Used to pass user options to the *this* to instruct it which data are desired.
1424 ***********************************************************************************************************/
1425
1426LUPI_HOST Distribution *parseGIDI( GIDI::Suite const &a_distribution, SetupInfo &a_setupInfo, Transporting::MC const &a_settings ) {
1427
1428 if( a_setupInfo.m_protare.projectileIntid( ) == PoPI::Intids::neutron ) {
1429 if( a_settings.wantRawTNSL_distributionSampling( ) ) {
1430 if( a_setupInfo.m_reaction->doubleDifferentialCrossSection( ).size( ) > 0 ) {
1431 GIDI::Form const *form = a_setupInfo.m_reaction->doubleDifferentialCrossSection( ).get<GIDI::Form>( 0 );
1432
1433 if( form->type( ) == GIDI::FormType::coherentElastic ) {
1436 return( new CoherentElasticTNSL( coherentElasticTNSL, a_setupInfo ) ); }
1437 else if( form->type( ) == GIDI::FormType::incoherentElastic ) {
1440 return( new IncoherentElasticTNSL( incoherentElasticTNSL, a_setupInfo ) );
1441 }
1442 }
1443 }
1444 }
1445
1446 std::string const *label = a_settings.styles( )->findLabelInLineage( a_distribution, a_setupInfo.m_distributionLabel );
1447 GIDI::Distributions::Distribution const &GIDI_distribution = *a_distribution.get<GIDI::Distributions::Distribution>( *label );
1448
1449 return parseGIDI2( GIDI_distribution, a_setupInfo, a_settings );
1450}
1451
1452/* *********************************************************************************************************//**
1453 * This function is used to convert the GIDI distribution *a_GIDI_distribution* into an MCGIDI distribution. It was split off of
1454 * **parseGIDI** to be able to handle referenced distribution by calling itself.
1455 *
1456 * @param a_distribution [in] The GIDI::Protare whose data is to be used to construct *this*.
1457 * @param a_setupInfo [in] Used internally when constructing a Protare to pass information to other constructors.
1458 * @param a_settings [in] Used to pass user options to the *this* to instruct it which data are desired.
1459 ***********************************************************************************************************/
1460
1461static LUPI_HOST Distribution *parseGIDI2( GIDI::Distributions::Distribution const &a_GIDI_distribution, SetupInfo &a_setupInfo,
1462 Transporting::MC const &a_settings ) {
1463
1464 Distribution *distribution = nullptr;
1465
1466 GIDI::FormType type = a_GIDI_distribution.type( );
1467
1468 switch( type ) {
1470 distribution = new AngularTwoBody( static_cast<GIDI::Distributions::AngularTwoBody const &>( a_GIDI_distribution ), a_setupInfo );
1471 break;
1473 distribution = new Uncorrelated( static_cast<GIDI::Distributions::Uncorrelated const &>( a_GIDI_distribution ), a_setupInfo );
1474 break;
1476 distribution = new KalbachMann( static_cast<GIDI::Distributions::KalbachMann const &>( a_GIDI_distribution ), a_setupInfo );
1477 break;
1479 distribution = new EnergyAngularMC( static_cast<GIDI::Distributions::EnergyAngularMC const &>( a_GIDI_distribution ), a_setupInfo );
1480 break;
1482 distribution = new AngularEnergyMC( static_cast<GIDI::Distributions::AngularEnergyMC const &>( a_GIDI_distribution ), a_setupInfo );
1483 break;
1485 distribution = new CoherentPhotoAtomicScattering( static_cast<GIDI::Distributions::CoherentPhotoAtomicScattering const &>( a_GIDI_distribution ), a_setupInfo );
1486 break;
1488 distribution = new IncoherentPhotoAtomicScattering( static_cast<GIDI::Distributions::IncoherentPhotoAtomicScattering const &>( a_GIDI_distribution ), a_setupInfo );
1489 break;
1491 distribution = new IncoherentBoundToFreePhotoAtomicScattering( static_cast<GIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering const &>( a_GIDI_distribution ), a_setupInfo );
1492 break;
1494 distribution = new Branching3d( static_cast<GIDI::Distributions::Branching3d const &>( a_GIDI_distribution ), a_setupInfo );
1495 break;
1497 distribution = new Unspecified( a_GIDI_distribution, a_setupInfo );
1498 break;
1500 GIDI::Distributions::Reference3d const *reference3d = static_cast<GIDI::Distributions::Reference3d const *>( &a_GIDI_distribution );
1501 GIDI::Distributions::Distribution const *linkedForm = static_cast<GIDI::Distributions::Distribution const *>( reference3d->findInAncestry( reference3d->href( ) ) );
1502 if( linkedForm == nullptr )
1503 throw std::runtime_error( "MCGIDI::Distributions::parseGIDI: could not find link '" + a_GIDI_distribution.toXLink( ) + "." );
1504 distribution = parseGIDI2( *linkedForm, a_setupInfo, a_settings ); }
1505 break;
1506 default :
1507 throw std::runtime_error( "MCGIDI::Distributions::parseGIDI: unsupported distribution: " + a_GIDI_distribution.toXLink( ) + "." );
1508 }
1509
1510 return( distribution );
1511}
1512
1513/* *********************************************************************************************************//**
1514 * @param a_distribution [in] The GIDI::Protare whose data is to be used to construct *this*.
1515 *
1516 * @return The type of the distribution or Distributions::Type::none if *a_distribution* is a *nullptr* pointer.
1517 ***********************************************************************************************************/
1518
1520
1521 if( a_distribution == nullptr ) return( Type::none );
1522 return( a_distribution->type( ) );
1523}
1524
1525}
1526
1527/* *********************************************************************************************************//**
1528 * This method serializes *this* for broadcasting as needed for MPI and GPUs. The method can count the number of required
1529 * bytes, pack *this* or unpack *this* depending on *a_mode*.
1530 *
1531 * @param a_buffer [in] The buffer to read or write data to depending on *a_mode*.A
1532 * @param a_mode [in] Specifies the action of this method.
1533 ***********************************************************************************************************/
1534
1536 Distributions::Distribution *a_distribution ) {
1537
1539 if( a_distribution != nullptr ) type = a_distribution->type( );
1540 int distributionType = distributionTypeToInt( type );
1541 DATA_MEMBER_INT( distributionType, a_buffer, a_mode );
1542 type = intToDistributionType( distributionType );
1543
1544 if( a_mode == LUPI::DataBuffer::Mode::Unpack ) {
1545 switch( type ) {
1547 a_distribution = nullptr;
1548 break;
1550 if (a_buffer.m_placement != nullptr) {
1551 a_distribution = new(a_buffer.m_placement) Distributions::Unspecified;
1552 a_buffer.incrementPlacement( sizeof( Distributions::Unspecified ) ); }
1553 else {
1554 a_distribution = new Distributions::Unspecified;
1555 }
1556 break;
1558 if (a_buffer.m_placement != nullptr) {
1559 a_distribution = new(a_buffer.m_placement) Distributions::AngularTwoBody;
1560 a_buffer.incrementPlacement( sizeof( Distributions::AngularTwoBody ) ); }
1561 else {
1562 a_distribution = new Distributions::AngularTwoBody;
1563 }
1564 break;
1566 if (a_buffer.m_placement != nullptr) {
1567 a_distribution = new(a_buffer.m_placement) Distributions::KalbachMann;
1568 a_buffer.incrementPlacement( sizeof( Distributions::KalbachMann ) ); }
1569 else {
1570 a_distribution = new Distributions::KalbachMann;
1571 }
1572 break;
1574 if (a_buffer.m_placement != nullptr) {
1575 a_distribution = new(a_buffer.m_placement) Distributions::Uncorrelated;
1576 a_buffer.incrementPlacement( sizeof( Distributions::Uncorrelated ) ); }
1577 else {
1578 a_distribution = new Distributions::Uncorrelated;
1579 }
1580 break;
1582 if (a_buffer.m_placement != nullptr) {
1583 a_distribution = new(a_buffer.m_placement) Distributions::Branching3d;
1584 a_buffer.incrementPlacement( sizeof( Distributions::Branching3d ) ); }
1585 else {
1586 a_distribution = new Distributions::Branching3d;
1587 }
1588 break;
1590 if (a_buffer.m_placement != nullptr) {
1591 a_distribution = new(a_buffer.m_placement) Distributions::EnergyAngularMC;
1592 a_buffer.incrementPlacement( sizeof( Distributions::EnergyAngularMC ) ); }
1593 else {
1594 a_distribution = new Distributions::EnergyAngularMC;
1595 }
1596 break;
1598 if (a_buffer.m_placement != nullptr) {
1599 a_distribution = new(a_buffer.m_placement) Distributions::AngularEnergyMC;
1600 a_buffer.incrementPlacement( sizeof( Distributions::AngularEnergyMC ) ); }
1601 else {
1602 a_distribution = new Distributions::AngularEnergyMC;
1603 }
1604 break;
1606 if( a_buffer.m_placement != nullptr ) {
1607 a_distribution = new(a_buffer.m_placement) Distributions::CoherentPhotoAtomicScattering;
1609 else {
1611 }
1612 break;
1614 if( a_buffer.m_placement != nullptr ) {
1615 a_distribution = new(a_buffer.m_placement) Distributions::IncoherentPhotoAtomicScattering;
1617 else {
1619 }
1620 break;
1622 if( a_buffer.m_placement != nullptr ) {
1625 else {
1627 }
1628 break;
1630 if( a_buffer.m_placement != nullptr ) {
1631 a_distribution = new(a_buffer.m_placement) Distributions::PairProductionGamma;
1633 else {
1634 a_distribution = new Distributions::PairProductionGamma;
1635 }
1636 break;
1638 if( a_buffer.m_placement != nullptr ) {
1639 a_distribution = new(a_buffer.m_placement) Distributions::CoherentElasticTNSL;
1641 else {
1642 a_distribution = new Distributions::CoherentElasticTNSL;
1643 }
1644 break;
1646 if( a_buffer.m_placement != nullptr ) {
1647 a_distribution = new(a_buffer.m_placement) Distributions::IncoherentElasticTNSL;
1649 else {
1650 a_distribution = new Distributions::IncoherentElasticTNSL;
1651 }
1652 break;
1654 if( a_buffer.m_placement != nullptr ) {
1657 else {
1659 }
1660 break;
1661 }
1662 }
1663
1664 if( a_mode == LUPI::DataBuffer::Mode::Memory ) {
1665 switch( type ) {
1667 break;
1669 a_buffer.incrementPlacement( sizeof( Distributions::Unspecified ) );
1670 break;
1673 break;
1675 a_buffer.incrementPlacement( sizeof( Distributions::KalbachMann ) );
1676 break;
1678 a_buffer.incrementPlacement( sizeof( Distributions::Uncorrelated ) );
1679 break;
1681 a_buffer.incrementPlacement( sizeof( Distributions::Branching3d ) );
1682 break;
1685 break;
1688 break;
1691 break;
1694 break;
1697 break;
1700 break;
1703 break;
1706 break;
1709 break;
1710 }
1711 }
1712
1713 switch( type ) {
1715 break;
1717 static_cast<Distributions::Unspecified *>( a_distribution )->serialize( a_buffer, a_mode );
1718 break;
1720 static_cast<Distributions::AngularTwoBody *>( a_distribution )->serialize( a_buffer, a_mode );
1721 break;
1723 static_cast<Distributions::KalbachMann *>( a_distribution )->serialize( a_buffer, a_mode );
1724 break;
1726 static_cast<Distributions::Uncorrelated *>( a_distribution )->serialize( a_buffer, a_mode );
1727 break;
1729 static_cast<Distributions::Branching3d *>( a_distribution )->serialize( a_buffer, a_mode );
1730 break;
1732 static_cast<Distributions::EnergyAngularMC *>( a_distribution )->serialize( a_buffer, a_mode );
1733 break;
1735 static_cast<Distributions::AngularEnergyMC *>( a_distribution )->serialize( a_buffer, a_mode );
1736 break;
1738 static_cast<Distributions::CoherentPhotoAtomicScattering *>( a_distribution )->serialize( a_buffer, a_mode );
1739 break;
1741 static_cast<Distributions::IncoherentPhotoAtomicScattering *>( a_distribution )->serialize( a_buffer, a_mode );
1742 break;
1744 static_cast<Distributions::IncoherentBoundToFreePhotoAtomicScattering *>( a_distribution )->serialize( a_buffer, a_mode );
1745 break;
1747 static_cast<Distributions::PairProductionGamma *>( a_distribution )->serialize( a_buffer, a_mode );
1748 break;
1750 static_cast<Distributions::CoherentElasticTNSL *>( a_distribution )->serialize( a_buffer, a_mode );
1751 break;
1753 static_cast<Distributions::IncoherentElasticTNSL *>( a_distribution )->serialize( a_buffer, a_mode );
1754 break;
1756 static_cast<Distributions::IncoherentPhotoAtomicScatteringElectron *>( a_distribution )->serialize( a_buffer, a_mode );
1757 break;
1758 }
1759
1760 return( a_distribution );
1761}
1762
1763}
G4double epsilon(G4double density, G4double temperature)
G4ThreadLocal T * G4GeomSplitter< T >::offset
#define GIDI_MonteCarlo_cdfStyleChars
Definition GIDI.hpp:245
#define DATA_MEMBER_VECTOR_DOUBLE(member, buf, mode)
#define DATA_MEMBER_DOUBLE(member, buf, mode)
#define DATA_MEMBER_INT( member, buf, mode)
#define LUPI_HOST_DEVICE
#define LUPI_HOST
#define PoPI_electronMass_MeV_c2
Definition PoPI.hpp:31
std::size_t size() const
Definition GIDI.hpp:900
std::vector< double > m_flattenedValues
Definition GIDI.hpp:898
std::string const & initialState() const
Definition GIDI.hpp:2492
FormType type() const
Definition GIDI.hpp:667
Axes const & axes() const
Definition GIDI.hpp:1012
Styles::Suite & styles()
Definition GIDI.hpp:4802
Component & doubleDifferentialCrossSection()
Definition GIDI.hpp:4297
std::string const * findLabelInLineage(GIDI::Suite const &a_suite, std::string const &a_label) const
T * get(std::size_t a_Index)
Definition GIDI.hpp:2642
std::size_t size() const
Definition GIDI.hpp:2591
Ancestry * findInAncestry(std::string const &a_href)
std::string toXLink() const
LUPI_HOST_DEVICE void incrementPlacement(std::size_t a_delta)
LUPI_HOST_DEVICE void serialize(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
LUPI_HOST_DEVICE Probabilities::ProbabilityBase2d_d1 * angular() const
LUPI_HOST_DEVICE void serialize(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
LUPI_HOST void setModelDBRC_data2(Sampling::Upscatter::ModelDBRC_data *a_modelDBRC_data)
LUPI_HOST_DEVICE Probabilities::ProbabilityBase2d_d1 * angular() const
LUPI_HOST_DEVICE void serialize(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
LUPI_HOST_DEVICE void serialize(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
LUPI_HOST_DEVICE double evaluate(double a_energyIn, double a_mu) const
LUPI_HOST_DEVICE void serialize(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
LUPI_HOST_DEVICE double evaluateFormFactor(double a_energyIn, double a_mu) const
LUPI_HOST_DEVICE void serialize(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
LUPI_HOST_DEVICE MCGIDI_VIRTUAL_FUNCTION ~Distribution() MCGIDI_TRUE_VIRTUAL
LUPI_HOST_DEVICE double targetMass() const
LUPI_HOST void setModelDBRC_data(Sampling::Upscatter::ModelDBRC_data *a_modelDBRC_data)
LUPI_HOST_DEVICE GIDI::Frame productFrame() const
LUPI_HOST_DEVICE double projectileMass() const
LUPI_HOST_DEVICE Probabilities::ProbabilityBase2d_d1 * energy() const
LUPI_HOST_DEVICE void serialize(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
LUPI_HOST_DEVICE double energyRatio(double a_energyIn, double a_mu) const
LUPI_HOST_DEVICE void serialize(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
LUPI_HOST_DEVICE double evaluateKleinNishina(double a_energyIn, double a_mu) const
LUPI_HOST_DEVICE double evaluateOccupationNumber(double a_X, double a_mu) const
LUPI_HOST_DEVICE void serialize(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
LUPI_HOST_DEVICE void serialize(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
LUPI_HOST_DEVICE double evaluateKleinNishina(double a_energyIn, double a_mu) const
LUPI_HOST_DEVICE void serialize(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
LUPI_HOST_DEVICE double energyRatio(double a_energyIn, double a_mu) const
LUPI_HOST_DEVICE double evaluateScatteringFactor(double a_X) const
LUPI_HOST_DEVICE Functions::Function2d * a() const
LUPI_HOST_DEVICE double evaluate(double E_in_lab, double E_out, double mu)
LUPI_HOST_DEVICE Functions::Function2d * r() const
LUPI_HOST_DEVICE void serialize(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
LUPI_HOST_DEVICE Probabilities::ProbabilityBase2d_d1 * f() const
LUPI_HOST_DEVICE void serialize(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
LUPI_HOST_DEVICE Probabilities::ProbabilityBase2d_d1 * angular() const
LUPI_HOST_DEVICE void serialize(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
LUPI_HOST_DEVICE Probabilities::ProbabilityBase2d * energy() const
LUPI_HOST_DEVICE void serialize(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
LUPI_HOST_DEVICE int projectileIntid() const
Definition MCGIDI.hpp:1517
GIDI::ProtareSingle const & m_GIDI_protare
Definition MCGIDI.hpp:255
std::string m_distributionLabel
Definition MCGIDI.hpp:269
std::map< std::string, int > m_stateNamesToIndices
Definition MCGIDI.hpp:277
ProtareSingle & m_protare
Definition MCGIDI.hpp:254
GIDI::Reaction const * m_reaction
Definition MCGIDI.hpp:272
LUPI_HOST GIDI::Styles::Suite const * styles() const
Definition MCGIDI.hpp:131
LUPI_HOST bool wantRawTNSL_distributionSampling() const
Definition MCGIDI.hpp:207
Definition GIDI.hpp:32
FormType
Definition GIDI.hpp:118
@ incoherentBoundToFreePhotonScattering
Definition GIDI.hpp:135
@ incoherentPhotonScattering
Definition GIDI.hpp:135
@ coherentPhotonScattering
Definition GIDI.hpp:135
@ incoherentElastic
Definition GIDI.hpp:136
Frame
Definition GIDI.hpp:146
@ centerOfMass
Definition GIDI.hpp:146
LUPI_HOST_DEVICE Type DistributionType(Distribution const *a_distribution)
LUPI_HOST Distribution * parseGIDI(GIDI::Suite const &a_distribution, SetupInfo &a_setupInfo, Transporting::MC const &a_settings)
LUPI_HOST Function1d_d1 * parseFunction1d_d1(Transporting::MC const &a_settings, GIDI::Suite const &a_suite)
Simple C++ string class, useful as replacement for std::string if this cannot be used,...
Definition MCGIDI.hpp:43
LUPI_HOST_DEVICE Distributions::Distribution * serializeDistribution(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode, Distributions::Distribution *a_distribution)
LUPI_HOST_DEVICE Probabilities::ProbabilityBase2d * serializeProbability2d(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode, Probabilities::ProbabilityBase2d *a_probability2d)
LUPI_HOST_DEVICE Functions::Function2d * serializeFunction2d(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode, Functions::Function2d *a_function2d)
LUPI_HOST_DEVICE Probabilities::ProbabilityBase2d_d1 * serializeProbability2d_d1(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode, Probabilities::ProbabilityBase2d_d1 *a_probability2d)
LUPI_HOST_DEVICE Probabilities::ProbabilityBase3d * serializeProbability3d(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode, Probabilities::ProbabilityBase3d *a_probability3d)
LUPI_HOST_DEVICE Interpolation GIDI2MCGIDI_interpolation(ptwXY_interpolation a_interpolation)
LUPI_HOST_DEVICE int binarySearchVector(double a_x, Vector< double > const &a_Xs, bool a_boundIndex=false)
Definition MCGIDI.hpp:318
LUPI_HOST_DEVICE int distributionTypeToInt(Distributions::Type a_type)
LUPI_HOST_DEVICE Functions::Function1d_d1 * serializeFunction1d_d1(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode, Functions::Function1d_d1 *a_function1d)
LUPI_HOST_DEVICE Distributions::Type intToDistributionType(int a_type)
ptwXY_interpolation ptwXY_stringToInterpolation(char const *interpolationString)
const axis_t axis_to_type< N >::axis
Definition pugixml.cc:9668
static int constexpr neutron
Definition PoPI.hpp:177