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