Geant4 Cross Reference |
1 /* 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 22 23 MCGIDI_outputChannel *outputChannel; 24 25 if( ( outputChannel = (MCGIDI_outputChanne 26 if( MCGIDI_outputChannel_initialize( smr, 27 return( outputChannel ); 28 } 29 /* 30 ********************************************** 31 */ 32 int MCGIDI_outputChannel_initialize( statusMes 33 34 memset( outputChannel, 0, sizeof( MCGIDI_o 35 return( 0 ); 36 } 37 /* 38 ********************************************** 39 */ 40 MCGIDI_outputChannel *MCGIDI_outputChannel_fre 41 42 MCGIDI_outputChannel_release( smr, outputC 43 smr_freeMemory( (void **) &outputChannel ) 44 return( NULL ); 45 } 46 /* 47 ********************************************** 48 */ 49 int MCGIDI_outputChannel_release( statusMessag 50 51 int i; 52 53 for( i = 0; i < outputChannel->numberOfPro 54 smr_freeMemory( (void **) &(outputChannel- 55 MCGIDI_outputChannel_initialize( smr, outp 56 57 return( 0 ); 58 } 59 /* 60 ********************************************** 61 */ 62 int MCGIDI_outputChannel_parseFromTOM( statusM 63 MCGIDI_reaction *reaction, MCGIDI_prod 64 65 int n{0}, delayedNeutronIndex{0}; 66 char const *genre{""}, *Q{""}; 67 xDataTOM_element *child{nullptr}; 68 69 MCGIDI_outputChannel_initialize( smr, outp 70 71 outputChannel->reaction = reaction; 72 outputChannel->parent = parent; 73 if( ( genre = xDataTOM_getAttributesValueI 74 if( ( parent != NULL ) && ( strcmp( genre, 75 smr_setReportError2( smr, smr_unknownI 76 goto err; 77 } 78 if( strcmp( genre, "twoBody" ) == 0 ) { 79 outputChannel->genre = MCGIDI_channelG 80 else if( strcmp( genre, "NBody" ) == 0 ) { 81 outputChannel->genre = MCGIDI_channelG 82 else if( strcmp( genre, "sumOfRemainingOut 83 outputChannel->genre = MCGIDI_channelG 84 else { 85 smr_setReportError2( smr, smr_unknownI 86 goto err; 87 } 88 if( ( Q = xDataTOM_getAttributesValueInEle 89 outputChannel->QIsFloat = !MCGIDI_misc_PQU 90 91 if( ( n = xDataTOM_numberOfElementsByName( 92 smr_setReportError2p( smr, smr_unknown 93 goto err; 94 } 95 if( ( outputChannel->products = (MCGIDI_pr 96 97 for( child = xDataTOME_getFirstElement( el 98 if( strcmp( child->name, "product" ) = 99 if( MCGIDI_product_parseFromTOM( s 100 &delayedNeutronIndex ) ) goto 101 outputChannel->numberOfProducts++; 102 else if( strcmp( child->name, "fission 103 continue; } 104 else { 105 printf( "outputChannel child not c 106 } 107 } 108 if( outputChannel->genre == MCGIDI_channel 109 double projectileMass_MeV, targetMass_ 110 111 projectileMass_MeV = MCGIDI_reaction_g 112 targetMass_MeV = MCGIDI_reaction_getTa 113 productMass_MeV = MCGIDI_product_getMa 114 residualMass_MeV = MCGIDI_product_getM 115 116 //TK 17-11-10 for v1.3 117 //A temporary fix for emission of gamm 118 // capture 119 if ( reaction->ENDF_MT == 102 && produ 120 //include/PoPs_data.h:#define e_Mas 121 residualMass_MeV += 5.4857990943e-4 122 } 123 124 MCGIDI_product_setTwoBodyMasses( smr, 125 } 126 127 return( 0 ); 128 129 err: 130 MCGIDI_outputChannel_release( smr, outputC 131 return( 1 ); 132 } 133 /* 134 ********************************************** 135 */ 136 int MCGIDI_outputChannel_numberOfProducts( MCG 137 138 return( outputChannel->numberOfProducts ); 139 } 140 /* 141 ********************************************** 142 */ 143 MCGIDI_product *MCGIDI_outputChannel_getProduc 144 145 if( ( i < 0 ) || ( i >= outputChannel->num 146 smr_setReportError2( smr, smr_unknownI 147 return( NULL ); 148 } 149 return( &(outputChannel->products[i]) ); 150 } 151 /* 152 ********************************************** 153 */ 154 int MCGIDI_outputChannel_getDomain( statusMess 155 156 if( outputChannel->reaction != NULL ) retu 157 return( MCGIDI_product_getDomain( smr, out 158 } 159 /* 160 ********************************************** 161 */ 162 MCGIDI_target_heated *MCGIDI_outputChannel_get 163 164 if( outputChannel->reaction != NULL ) retu 165 return( MCGIDI_product_getTargetHeated( sm 166 } 167 /* 168 ********************************************** 169 */ 170 double MCGIDI_outputChannel_getProjectileMass_ 171 172 if( outputChannel->reaction != NULL ) retu 173 return( MCGIDI_product_getProjectileMass_M 174 } 175 /* 176 ********************************************** 177 */ 178 double MCGIDI_outputChannel_getTargetMass_MeV( 179 180 if( outputChannel->reaction != NULL ) retu 181 return( MCGIDI_product_getTargetMass_MeV( 182 } 183 /* 184 ********************************************** 185 */ 186 double MCGIDI_outputChannel_getQ_MeV( statusMe 187 188 return( outputChannel->Q ); 189 } 190 /* 191 ********************************************** 192 */ 193 double MCGIDI_outputChannel_getFinalQ( statusM 194 195 int iProduct; 196 double Q = outputChannel->Q; 197 MCGIDI_product *product; 198 199 for( iProduct = 0; iProduct < outputChanne 200 product = &(outputChannel->products[iP 201 if( product->decayChannel.genre != MCG 202 if( !smr_isOk( smr ) ) break; 203 } 204 return( Q ); 205 } 206 /* 207 ********************************************** 208 */ 209 int MCGIDI_outputChannel_sampleProductsAtE(sta 210 MCG 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 293 case MCGIDI_channelGenre_sumOfRemain 294 // Get mass of final state particl 295 // masses[0] and masses[1] are inc 296 masses[2] = MCGIDI_product_getMass 297 switch( distribution->type ) { 298 case MCGIDI_distributionType_uncor 299 MCGIDI_uncorrelated_sampleDistri 300 break; 301 case MCGIDI_distributionType_energ 302 MCGIDI_energyAngular_sampleDistr 303 break; 304 case MCGIDI_distributionType_Kalba 305 MCGIDI_KalbachMann_sampleEp( smr 306 break; 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 316 case MCGIDI_channelGenre_undefined_e 317 printf( "Channel is undefined\n" ) 318 break; 319 320 case MCGIDI_channelGenre_twoBodyDeca 321 printf( "Channel is twoBodyDecay\n 322 break; 323 324 case MCGIDI_channelGenre_uncorrelate 325 printf( "Channel is uncorrelatedDe 326 break; 327 328 default : 329 printf( "Unsupported channel genre 330 break; 331 } 332 333 if (!smr_isOk(smr) ) return( -1 ); 334 if (!secondTwoBody) { 335 if (decaySamplingInfo->frame == xD 336 if (MCGIDI_kinetics_COM2Lab( smr 337 } 338 339 // Assign kinematics to final stat 340 productData[0].kineticEnergy = dec 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 } 360 361 #if defined __cplusplus 362 } 363 #endif 364 365