| Geant4 Cross Reference |
1 /*
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5 #include <string.h>
6 #define _USE_MATH_DEFINES
7 #include <cmath>
8
9 #include "MCGIDI.h"
10 #include "MCGIDI_misc.h"
11
12 #if defined __cplusplus
13 #include "G4Log.hh"
14 namespace GIDI {
15 using namespace GIDI;
16 #endif
17
18 /*
19 ************************************************************
20 */
21 MCGIDI_outputChannel *MCGIDI_outputChannel_new( statusMessageReporting *smr ) {
22
23 MCGIDI_outputChannel *outputChannel;
24
25 if( ( outputChannel = (MCGIDI_outputChannel *) smr_malloc2( smr, sizeof( MCGIDI_outputChannel ), 0, "outputChannel" ) ) == NULL ) return( NULL );
26 if( MCGIDI_outputChannel_initialize( smr, outputChannel ) ) outputChannel = MCGIDI_outputChannel_free( smr, outputChannel );
27 return( outputChannel );
28 }
29 /*
30 ************************************************************
31 */
32 int MCGIDI_outputChannel_initialize( statusMessageReporting * /*smr*/, MCGIDI_outputChannel *outputChannel ) {
33
34 memset( outputChannel, 0, sizeof( MCGIDI_outputChannel ) );
35 return( 0 );
36 }
37 /*
38 ************************************************************
39 */
40 MCGIDI_outputChannel *MCGIDI_outputChannel_free( statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel ) {
41
42 MCGIDI_outputChannel_release( smr, outputChannel );
43 smr_freeMemory( (void **) &outputChannel );
44 return( NULL );
45 }
46 /*
47 ************************************************************
48 */
49 int MCGIDI_outputChannel_release( statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel ) {
50
51 int i;
52
53 for( i = 0; i < outputChannel->numberOfProducts; i++ ) MCGIDI_product_release( smr, &(outputChannel->products[i]) );
54 smr_freeMemory( (void **) &(outputChannel->products) );
55 MCGIDI_outputChannel_initialize( smr, outputChannel );
56
57 return( 0 );
58 }
59 /*
60 ************************************************************
61 */
62 int MCGIDI_outputChannel_parseFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_POPs *pops, MCGIDI_outputChannel *outputChannel,
63 MCGIDI_reaction *reaction, MCGIDI_product *parent ) {
64
65 int n{0}, delayedNeutronIndex{0};
66 char const *genre{""}, *Q{""};
67 xDataTOM_element *child{nullptr};
68
69 MCGIDI_outputChannel_initialize( smr, outputChannel );
70
71 outputChannel->reaction = reaction;
72 outputChannel->parent = parent;
73 if( ( genre = xDataTOM_getAttributesValueInElement( element, "genre" ) ) == NULL ) goto err;
74 if( ( parent != NULL ) && ( strcmp( genre, "NBody" ) ) ) {
75 smr_setReportError2( smr, smr_unknownID, 1, "decay channel's genre can only be 'uncorreclated' (a.k.a. 'NBody') and not '%s'", genre );
76 goto err;
77 }
78 if( strcmp( genre, "twoBody" ) == 0 ) {
79 outputChannel->genre = MCGIDI_channelGenre_twoBody_e; }
80 else if( strcmp( genre, "NBody" ) == 0 ) {
81 outputChannel->genre = MCGIDI_channelGenre_uncorrelated_e; }
82 else if( strcmp( genre, "sumOfRemainingOutputChannels" ) == 0 ) {
83 outputChannel->genre = MCGIDI_channelGenre_sumOfRemaining_e; }
84 else {
85 smr_setReportError2( smr, smr_unknownID, 1, "unsupported genre = '%s'", genre );
86 goto err;
87 }
88 if( ( Q = xDataTOM_getAttributesValueInElement( element, "Q" ) ) == NULL ) goto err;
89 outputChannel->QIsFloat = !MCGIDI_misc_PQUStringToDoubleInUnitOf( smr, Q, "MeV", &(outputChannel->Q) );
90
91 if( ( n = xDataTOM_numberOfElementsByName( smr, element, "product" ) ) == 0 ) {
92 smr_setReportError2p( smr, smr_unknownID, 1, "outputChannel does not have any products" );
93 goto err;
94 }
95 if( ( outputChannel->products = (MCGIDI_product *) smr_malloc2( smr, n * sizeof( MCGIDI_product ), 0, "outputChannel->products" ) ) == NULL ) goto err;
96
97 for( child = xDataTOME_getFirstElement( element ); child != NULL; child = xDataTOME_getNextElement( child ) ) {
98 if( strcmp( child->name, "product" ) == 0 ) {
99 if( MCGIDI_product_parseFromTOM( smr, child, outputChannel, pops, &(outputChannel->products[outputChannel->numberOfProducts]),
100 &delayedNeutronIndex ) ) goto err;
101 outputChannel->numberOfProducts++; }
102 else if( strcmp( child->name, "fissionEnergyReleased" ) == 0 ) { /* ????????? Need to support. */
103 continue; }
104 else {
105 printf( "outputChannel child not currently supported = %s\n", child->name );
106 }
107 }
108 if( outputChannel->genre == MCGIDI_channelGenre_twoBody_e ) {
109 double projectileMass_MeV, targetMass_MeV, productMass_MeV, residualMass_MeV;
110
111 projectileMass_MeV = MCGIDI_reaction_getProjectileMass_MeV( smr, reaction );
112 targetMass_MeV = MCGIDI_reaction_getTargetMass_MeV( smr, reaction );
113 productMass_MeV = MCGIDI_product_getMass_MeV( smr, &(outputChannel->products[0]) );
114 residualMass_MeV = MCGIDI_product_getMass_MeV( smr, &(outputChannel->products[1]) );
115
116 //TK 17-11-10 for v1.3
117 //A temporary fix for emission of gamma(2.2MeV) from n captured by H
118 // capture gamma D
119 if ( reaction->ENDF_MT == 102 && productMass_MeV == 0 && ( outputChannel->products[1].pop->A == 2 && outputChannel->products[1].pop->Z == 1 ) ) {
120 //include/PoPs_data.h:#define e_Mass 5.4857990943e-4 /* electron mass in AMU */
121 residualMass_MeV += 5.4857990943e-4*MCGIDI_AMU2MeV;
122 }
123
124 MCGIDI_product_setTwoBodyMasses( smr, &(outputChannel->products[0]), projectileMass_MeV, targetMass_MeV, productMass_MeV, residualMass_MeV );
125 }
126
127 return( 0 );
128
129 err:
130 MCGIDI_outputChannel_release( smr, outputChannel );
131 return( 1 );
132 }
133 /*
134 ************************************************************
135 */
136 int MCGIDI_outputChannel_numberOfProducts( MCGIDI_outputChannel *outputChannel ) {
137
138 return( outputChannel->numberOfProducts );
139 }
140 /*
141 ************************************************************
142 */
143 MCGIDI_product *MCGIDI_outputChannel_getProductAtIndex( statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel, int i ) {
144
145 if( ( i < 0 ) || ( i >= outputChannel->numberOfProducts ) ) {
146 smr_setReportError2( smr, smr_unknownID, 1, "bad product index = %d: outputChannel as only %d products", i, outputChannel->numberOfProducts );
147 return( NULL );
148 }
149 return( &(outputChannel->products[i]) );
150 }
151 /*
152 ************************************************************
153 */
154 int MCGIDI_outputChannel_getDomain( statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel, double *EMin, double *EMax ) {
155
156 if( outputChannel->reaction != NULL ) return( MCGIDI_reaction_getDomain( smr, outputChannel->reaction, EMin, EMax ) );
157 return( MCGIDI_product_getDomain( smr, outputChannel->parent, EMin, EMax ) );
158 }
159 /*
160 ************************************************************
161 */
162 MCGIDI_target_heated *MCGIDI_outputChannel_getTargetHeated( statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel ) {
163
164 if( outputChannel->reaction != NULL ) return( MCGIDI_reaction_getTargetHeated( smr, outputChannel->reaction ) );
165 return( MCGIDI_product_getTargetHeated( smr, outputChannel->parent ) );
166 }
167 /*
168 ************************************************************
169 */
170 double MCGIDI_outputChannel_getProjectileMass_MeV( statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel ) {
171
172 if( outputChannel->reaction != NULL ) return( MCGIDI_reaction_getProjectileMass_MeV( smr, outputChannel->reaction ) );
173 return( MCGIDI_product_getProjectileMass_MeV( smr, outputChannel->parent ) );
174 }
175 /*
176 ************************************************************
177 */
178 double MCGIDI_outputChannel_getTargetMass_MeV( statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel ) {
179
180 if( outputChannel->reaction != NULL ) return( MCGIDI_reaction_getTargetMass_MeV( smr, outputChannel->reaction ) );
181 return( MCGIDI_product_getTargetMass_MeV( smr, outputChannel->parent ) );
182 }
183 /*
184 ************************************************************
185 */
186 double MCGIDI_outputChannel_getQ_MeV( statusMessageReporting * /*smr*/, MCGIDI_outputChannel *outputChannel, double /*e_in*/ ) {
187
188 return( outputChannel->Q );
189 }
190 /*
191 ************************************************************
192 */
193 double MCGIDI_outputChannel_getFinalQ( statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel, double e_in ) {
194
195 int iProduct;
196 double Q = outputChannel->Q;
197 MCGIDI_product *product;
198
199 for( iProduct = 0; iProduct < outputChannel->numberOfProducts; iProduct++ ) {
200 product = &(outputChannel->products[iProduct]);
201 if( product->decayChannel.genre != MCGIDI_channelGenre_undefined_e ) Q += MCGIDI_outputChannel_getFinalQ( smr, &(product->decayChannel), e_in );
202 if( !smr_isOk( smr ) ) break;
203 }
204 return( Q );
205 }
206 /*
207 ************************************************************
208 */
209 int MCGIDI_outputChannel_sampleProductsAtE(statusMessageReporting* smr,
210 MCGIDI_outputChannel* outputChannel,
211 MCGIDI_quantitiesLookupModes &modes,
212 MCGIDI_decaySamplingInfo* decaySamplingInfo,
213 MCGIDI_sampledProductsDatas* productDatas,
214 double *masses_ )
215 {
216 int i1;
217 int multiplicity(0);
218 int secondTwoBody = 0, isDecayChannel = ( outputChannel->reaction == NULL );
219 double e_in = modes.getProjectileEnergy( );
220 MCGIDI_product *product;
221 double phi, p, masses[3];
222 MCGIDI_distribution *distribution;
223 MCGIDI_sampledProductsData productData[2];
224
225 if (isDecayChannel) {
226 masses[0] = masses_[0]; /* More work may be needed here. */
227 masses[1] = masses_[1];
228 } else {
229 masses[0] = MCGIDI_reaction_getProjectileMass_MeV( smr, outputChannel->reaction );
230 masses[1] = MCGIDI_reaction_getTargetMass_MeV( smr, outputChannel->reaction );
231 }
232
233 // Loop over all possible final state particles reachable from initial state
234 // List of these particles (products) was read in from GIDI
235 // Note: all particles satifying the sampling criteria are included in the
236 // final state, regardless of charge, energy or baryon number conservation
237
238 for (i1 = 0; i1 < outputChannel->numberOfProducts; i1++) {
239 product = &(outputChannel->products[i1]);
240 if (product->decayChannel.genre != MCGIDI_channelGenre_undefined_e ) {
241 if( MCGIDI_outputChannel_sampleProductsAtE(smr, &(product->decayChannel),
242 modes, decaySamplingInfo,
243 productDatas, masses ) < 0 ) return( -1 );
244 } else {
245 distribution = &(product->distribution);
246 if( distribution->type == MCGIDI_distributionType_none_e ) continue;
247
248 if (!secondTwoBody) {
249 // Sample multiplicity of final state particle at kinetic energy of projectile
250 // The multiplicity stored in GIDI is a real number whose fractional part is
251 // compared to a random number to decide what integer value is returned
252 if ((multiplicity = product->multiplicity) == 0) multiplicity =
253 MCGIDI_product_sampleMultiplicity(smr, product, e_in,
254 decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
255 while (multiplicity > 0) {
256
257 multiplicity--;
258 decaySamplingInfo->pop = product->pop;
259 decaySamplingInfo->mu = 0;
260 decaySamplingInfo->Ep = 0;
261 productData[0].isVelocity = decaySamplingInfo->isVelocity;
262 productData[0].pop = product->pop;
263 productData[0].delayedNeutronIndex = product->delayedNeutronIndex;
264 productData[0].delayedNeutronRate = product->delayedNeutronRate;
265 productData[0].birthTimeSec = 0;
266 if (product->delayedNeutronRate > 0) {
267 productData[0].birthTimeSec =
268 -G4Log( decaySamplingInfo->rng( decaySamplingInfo->rngState ) ) / product->delayedNeutronRate;
269 }
270
271 switch( outputChannel->genre ) {
272
273 case MCGIDI_channelGenre_twoBody_e :
274 secondTwoBody = 1;
275 MCGIDI_angular_sampleMu( smr, distribution->angular, modes, decaySamplingInfo );
276 if (smr_isOk(smr) ) {
277 phi = 2. * M_PI * decaySamplingInfo->rng( decaySamplingInfo->rngState );
278 MCGIDI_kinetics_2BodyReaction( smr, distribution->angular, e_in, decaySamplingInfo->mu, phi, productData );
279 if (!smr_isOk(smr) ) return( -1 );
280 productData[1].pop = product[1].pop;
281 productData[1].delayedNeutronIndex = product[1].delayedNeutronIndex;
282 productData[1].delayedNeutronRate = product->delayedNeutronRate;
283 productData[1].birthTimeSec = 0;
284 MCGIDI_sampledProducts_addProduct( smr, productDatas, productData );
285 if( !smr_isOk( smr ) ) return( -1 );
286 MCGIDI_sampledProducts_addProduct( smr, productDatas, &(productData[1]) );
287 if( !smr_isOk( smr ) ) return( -1 );
288 }
289 break;
290
291 case MCGIDI_channelGenre_uncorrelated_e :
292
293 case MCGIDI_channelGenre_sumOfRemaining_e :
294 // Get mass of final state particle, then get its distribution
295 // masses[0] and masses[1] are incident and target masses
296 masses[2] = MCGIDI_product_getMass_MeV( smr, product );
297 switch( distribution->type ) {
298 case MCGIDI_distributionType_uncorrelated_e :
299 MCGIDI_uncorrelated_sampleDistribution( smr, distribution, modes, decaySamplingInfo );
300 break;
301 case MCGIDI_distributionType_energyAngular_e :
302 MCGIDI_energyAngular_sampleDistribution( smr, distribution, modes, decaySamplingInfo );
303 break;
304 case MCGIDI_distributionType_KalbachMann_e :
305 MCGIDI_KalbachMann_sampleEp( smr, distribution->KalbachMann, modes, decaySamplingInfo );
306 break;
307 case MCGIDI_distributionType_angularEnergy_e :
308 MCGIDI_angularEnergy_sampleDistribution( smr, distribution->angularEnergy, modes, decaySamplingInfo );
309 break;
310 default :
311 printf( "Unknown spectral data form product name = %s, channel genre = %d\n", product->pop->name, outputChannel->genre );
312 break;
313 }
314 break;
315
316 case MCGIDI_channelGenre_undefined_e :
317 printf( "Channel is undefined\n" );
318 break;
319
320 case MCGIDI_channelGenre_twoBodyDecay_e :
321 printf( "Channel is twoBodyDecay\n" );
322 break;
323
324 case MCGIDI_channelGenre_uncorrelatedDecay_e :
325 printf( "Channel is uncorrelatedDecay\n" );
326 break;
327
328 default :
329 printf( "Unsupported channel genre = %d\n", outputChannel->genre );
330 break;
331 }
332
333 if (!smr_isOk(smr) ) return( -1 );
334 if (!secondTwoBody) {
335 if (decaySamplingInfo->frame == xDataTOM_frame_centerOfMass) {
336 if (MCGIDI_kinetics_COM2Lab( smr, modes, decaySamplingInfo, masses) != 0 ) return( -1 );
337 }
338
339 // Assign kinematics to final state product
340 productData[0].kineticEnergy = decaySamplingInfo->Ep;
341 p = std::sqrt( decaySamplingInfo->Ep * ( decaySamplingInfo->Ep + 2. * product->pop->mass_MeV ) );
342 if (productData[0].isVelocity) p *= MCGIDI_speedOfLight_cm_sec / std::sqrt( p * p + product->pop->mass_MeV * product->pop->mass_MeV );
343 productData[0].pz_vz = p * decaySamplingInfo->mu;
344 p = std::sqrt( 1. - decaySamplingInfo->mu * decaySamplingInfo->mu ) * p;
345 phi = 2. * M_PI * decaySamplingInfo->rng( decaySamplingInfo->rngState );
346 productData[0].px_vx = p * std::sin( phi );
347 productData[0].py_vy = p * std::cos( phi );
348 MCGIDI_sampledProducts_addProduct( smr, productDatas, productData );
349 if (!smr_isOk(smr) ) return( -1 );
350 }
351 } // while multiplicity
352
353 } // if !secondTwoBody
354 } // if decay channel genre
355
356 } // loop over possible final state products
357 return( productDatas->numberOfProducts );
358
359 }
360
361 #if defined __cplusplus
362 }
363 #endif
364
365