9 #include "MCGIDI_misc.h"
11 #if defined __cplusplus
22 MCGIDI_outputChannel *outputChannel;
24 if( ( outputChannel = (MCGIDI_outputChannel *) smr_malloc2( smr,
sizeof( MCGIDI_outputChannel ), 0,
"outputChannel" ) ) == NULL )
return( NULL );
26 return( outputChannel );
33 memset( outputChannel, 0,
sizeof( MCGIDI_outputChannel ) );
52 for( i = 0; i < outputChannel->numberOfProducts; i++ )
MCGIDI_product_release( smr, &(outputChannel->products[i]) );
62 MCGIDI_reaction *reaction, MCGIDI_product *parent ) {
64 int n, delayedNeutronIndex = 0;
65 char const *genre, *
Q;
66 xDataTOM_element *child;
70 outputChannel->reaction = reaction;
71 outputChannel->parent = parent;
73 if( ( parent != NULL ) && ( strcmp( genre,
"NBody" ) ) ) {
74 smr_setReportError2( smr, smr_unknownID, 1,
"decay channel's genre can only be 'uncorreclated' (a.k.a. 'NBody') and not '%s'", genre );
77 if( strcmp( genre,
"twoBody" ) == 0 ) {
78 outputChannel->genre = MCGIDI_channelGenre_twoBody_e; }
79 else if( strcmp( genre,
"NBody" ) == 0 ) {
80 outputChannel->genre = MCGIDI_channelGenre_uncorrelated_e; }
81 else if( strcmp( genre,
"sumOfRemainingOutputChannels" ) == 0 ) {
82 outputChannel->genre = MCGIDI_channelGenre_sumOfRemaining_e; }
84 smr_setReportError2( smr, smr_unknownID, 1,
"unsupported genre = '%s'", genre );
91 smr_setReportError2p( smr, smr_unknownID, 1,
"outputChannel does not have any products" );
94 if( ( outputChannel->products = (MCGIDI_product *) smr_malloc2( smr, n *
sizeof( MCGIDI_product ), 0,
"outputChannel->products" ) ) == NULL )
goto err;
97 if( strcmp( child->name,
"product" ) == 0 ) {
99 &delayedNeutronIndex ) )
goto err;
100 outputChannel->numberOfProducts++; }
101 else if( strcmp( child->name,
"fissionEnergyReleased" ) == 0 ) {
104 printf(
"outputChannel child not currently supported = %s\n", child->name );
107 if( outputChannel->genre == MCGIDI_channelGenre_twoBody_e ) {
108 double projectileMass_MeV, targetMass_MeV, productMass_MeV, residualMass_MeV;
128 return( outputChannel->numberOfProducts );
135 if( ( i < 0 ) || ( i >= outputChannel->numberOfProducts ) ) {
136 smr_setReportError2( smr, smr_unknownID, 1,
"bad product index = %d: outputChannel as only %d products", i, outputChannel->numberOfProducts );
139 return( &(outputChannel->products[i]) );
178 return( outputChannel->Q );
186 double Q = outputChannel->Q;
187 MCGIDI_product *product;
189 for( iProduct = 0; iProduct < outputChannel->numberOfProducts; iProduct++ ) {
190 product = &(outputChannel->products[iProduct]);
200 MCGIDI_decaySamplingInfo *decaySamplingInfo, MCGIDI_sampledProductsDatas *productDatas,
double *masses_ ) {
202 int i1, multiplicity, secondTwoBody = 0, isDecayChannel = ( outputChannel->reaction == NULL );
203 double e_in = modes.getProjectileEnergy( );
204 MCGIDI_product *product;
205 double phi, p, masses[3];
206 MCGIDI_distribution *distribution;
207 MCGIDI_sampledProductsData productData[2];
209 if( isDecayChannel ) {
210 masses[0] = masses_[0];
211 masses[1] = masses_[1]; }
217 for( i1 = 0; i1 < outputChannel->numberOfProducts; i1++ ) {
218 product = &(outputChannel->products[i1]);
219 if( product->decayChannel.genre != MCGIDI_channelGenre_undefined_e ) {
222 distribution = &(product->distribution);
223 if( distribution->type == MCGIDI_distributionType_none_e )
continue;
224 if( !secondTwoBody ) {
226 decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
227 while( multiplicity > 0 ) {
230 decaySamplingInfo->pop = product->pop;
231 decaySamplingInfo->mu = 0;
232 decaySamplingInfo->Ep = 0;
233 productData[0].isVelocity = decaySamplingInfo->isVelocity;
234 productData[0].pop = product->pop;
235 productData[0].delayedNeutronIndex = product->delayedNeutronIndex;
236 productData[0].delayedNeutronRate = product->delayedNeutronRate;
237 productData[0].birthTimeSec = 0;
238 if( product->delayedNeutronRate > 0 ) {
239 productData[0].birthTimeSec = -
G4Log( decaySamplingInfo->rng( decaySamplingInfo->rngState ) ) / product->delayedNeutronRate;
242 switch( outputChannel->genre ) {
243 case MCGIDI_channelGenre_twoBody_e :
247 phi = 2. * M_PI * decaySamplingInfo->rng( decaySamplingInfo->rngState );
249 if( !
smr_isOk( smr ) )
return( -1 );
250 productData[1].pop = product[1].pop;
251 productData[1].delayedNeutronIndex = product[1].delayedNeutronIndex;
252 productData[1].delayedNeutronRate = product->delayedNeutronRate;
253 productData[1].birthTimeSec = 0;
255 if( !
smr_isOk( smr ) )
return( -1 );
257 if( !
smr_isOk( smr ) )
return( -1 );
260 case MCGIDI_channelGenre_uncorrelated_e :
261 case MCGIDI_channelGenre_sumOfRemaining_e :
263 switch( distribution->type ) {
264 case MCGIDI_distributionType_uncorrelated_e :
267 case MCGIDI_distributionType_energyAngular_e :
270 case MCGIDI_distributionType_KalbachMann_e :
273 case MCGIDI_distributionType_angularEnergy_e :
277 printf(
"Unknown spectral data form product name = %s, channel genre = %d\n", product->pop->name, outputChannel->genre );
281 case MCGIDI_channelGenre_undefined_e :
282 printf(
"Channel is undefined\n" );
283 case MCGIDI_channelGenre_twoBodyDecay_e :
284 printf(
"Channel is twoBodyDecay\n" );
285 case MCGIDI_channelGenre_uncorrelatedDecay_e :
286 printf(
"Channel is uncorrelatedDecay\n" );
288 printf(
"Unsupported channel genre = %d\n", outputChannel->genre );
290 if( !
smr_isOk( smr ) )
return( -1 );
291 if( !secondTwoBody ) {
292 if( decaySamplingInfo->frame == xDataTOM_frame_centerOfMass ) {
295 productData[0].kineticEnergy = decaySamplingInfo->Ep;
296 p = std::sqrt( decaySamplingInfo->Ep * ( decaySamplingInfo->Ep + 2. * product->pop->mass_MeV ) );
297 if( productData[0].isVelocity ) p *= MCGIDI_speedOfLight_cm_sec / std::sqrt( p * p + product->pop->mass_MeV * product->pop->mass_MeV );
298 productData[0].pz_vz = p * decaySamplingInfo->mu;
299 p = std::sqrt( 1. - decaySamplingInfo->mu * decaySamplingInfo->mu ) * p;
300 phi = 2. * M_PI * decaySamplingInfo->rng( decaySamplingInfo->rngState );
301 productData[0].px_vx = p * std::sin( phi );
302 productData[0].py_vy = p * std::cos( phi );
304 if( !
smr_isOk( smr ) )
return( -1 );
310 return( productDatas->numberOfProducts );
313 #if defined __cplusplus
int MCGIDI_product_parseFromTOM(statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_outputChannel *outputChannel, MCGIDI_POPs *pops, MCGIDI_product *product, int *delayedNeutronIndex)
int MCGIDI_kinetics_COM2Lab(statusMessageReporting *smr, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo, double masses[3])
int MCGIDI_angular_sampleMu(statusMessageReporting *smr, MCGIDI_angular *angular, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)
int MCGIDI_outputChannel_release(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel)
int MCGIDI_product_sampleMultiplicity(statusMessageReporting *, MCGIDI_product *product, double e_in, double r)
xDataTOM_element * xDataTOME_getFirstElement(xDataTOM_element *element)
int MCGIDI_reaction_getDomain(statusMessageReporting *, MCGIDI_reaction *reaction, double *EMin, double *EMax)
double MCGIDI_outputChannel_getQ_MeV(statusMessageReporting *, MCGIDI_outputChannel *outputChannel, double)
double MCGIDI_reaction_getTargetMass_MeV(statusMessageReporting *smr, MCGIDI_reaction *reaction)
int MCGIDI_KalbachMann_sampleEp(statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)
double MCGIDI_product_getTargetMass_MeV(statusMessageReporting *smr, MCGIDI_product *product)
int MCGIDI_kinetics_2BodyReaction(statusMessageReporting *smr, MCGIDI_angular *angular, double K, double mu, double phi, MCGIDI_sampledProductsData *outgoingData)
int xDataTOM_numberOfElementsByName(statusMessageReporting *, xDataTOM_element *element, char const *name)
double MCGIDI_product_getMass_MeV(statusMessageReporting *, MCGIDI_product *product)
double MCGIDI_outputChannel_getTargetMass_MeV(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel)
double MCGIDI_outputChannel_getFinalQ(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel, double e_in)
int MCGIDI_angularEnergy_sampleDistribution(statusMessageReporting *smr, MCGIDI_angularEnergy *angularEnergy, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)
int MCGIDI_product_release(statusMessageReporting *smr, MCGIDI_product *product)
int smr_isOk(statusMessageReporting *smr)
int MCGIDI_outputChannel_parseFromTOM(statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_POPs *pops, MCGIDI_outputChannel *outputChannel, MCGIDI_reaction *reaction, MCGIDI_product *parent)
MCGIDI_target_heated * MCGIDI_reaction_getTargetHeated(statusMessageReporting *, MCGIDI_reaction *reaction)
void * smr_freeMemory(void **p)
MCGIDI_product * MCGIDI_outputChannel_getProductAtIndex(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel, int i)
MCGIDI_target_heated * MCGIDI_outputChannel_getTargetHeated(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel)
int MCGIDI_outputChannel_initialize(statusMessageReporting *, MCGIDI_outputChannel *outputChannel)
int MCGIDI_sampledProducts_addProduct(statusMessageReporting *smr, MCGIDI_sampledProductsDatas *sampledProductsDatas, MCGIDI_sampledProductsData *sampledProductsData)
double MCGIDI_product_getProjectileMass_MeV(statusMessageReporting *smr, MCGIDI_product *product)
int MCGIDI_misc_PQUStringToDoubleInUnitOf(statusMessageReporting *smr, char const *str, char const *toUnit, double *value)
MCGIDI_target_heated * MCGIDI_product_getTargetHeated(statusMessageReporting *smr, MCGIDI_product *product)
char const * xDataTOM_getAttributesValueInElement(xDataTOM_element *element, char const *name)
int MCGIDI_outputChannel_sampleProductsAtE(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo, MCGIDI_sampledProductsDatas *productDatas, double *masses_)
G4double G4Log(G4double x)
double MCGIDI_reaction_getProjectileMass_MeV(statusMessageReporting *smr, MCGIDI_reaction *reaction)
xDataTOM_element * xDataTOME_getNextElement(xDataTOM_element *element)
MCGIDI_outputChannel * MCGIDI_outputChannel_free(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel)
int MCGIDI_outputChannel_getDomain(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel, double *EMin, double *EMax)
int MCGIDI_outputChannel_numberOfProducts(MCGIDI_outputChannel *outputChannel)
MCGIDI_outputChannel * MCGIDI_outputChannel_new(statusMessageReporting *smr)
int MCGIDI_product_getDomain(statusMessageReporting *smr, MCGIDI_product *product, double *EMin, double *EMax)
int MCGIDI_product_setTwoBodyMasses(statusMessageReporting *smr, MCGIDI_product *product, double projectileMass_MeV, double targetMass_MeV, double productMass_MeV, double residualMass_MeV)
double MCGIDI_outputChannel_getProjectileMass_MeV(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel)
int MCGIDI_energyAngular_sampleDistribution(statusMessageReporting *smr, MCGIDI_distribution *distribution, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)