Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 /* 26 /* 27 # <<BEGIN-copyright>> 27 # <<BEGIN-copyright>> 28 # <<END-copyright>> 28 # <<END-copyright>> 29 */ 29 */ 30 #include <iostream> 30 #include <iostream> 31 #include <stdlib.h> 31 #include <stdlib.h> 32 #include <cmath> 32 #include <cmath> 33 33 34 #include "G4GIDI_target.hh" 34 #include "G4GIDI_target.hh" 35 #include "G4GIDI_mass.hh" 35 #include "G4GIDI_mass.hh" 36 #include "G4GIDI_Misc.hh" 36 #include "G4GIDI_Misc.hh" 37 37 38 using namespace std; 38 using namespace std; 39 using namespace GIDI; 39 using namespace GIDI; 40 40 41 /* 41 /* 42 ********************************************** 42 *************************************************************** 43 */ 43 */ 44 G4GIDI_target::G4GIDI_target( char const *file 44 G4GIDI_target::G4GIDI_target( char const *fileName ) { 45 45 46 init( fileName ); 46 init( fileName ); 47 } 47 } 48 /* 48 /* 49 ********************************************** 49 *************************************************************** 50 */ 50 */ 51 G4GIDI_target::G4GIDI_target( string const &fi 51 G4GIDI_target::G4GIDI_target( string const &fileName ) { 52 52 53 init( fileName.c_str( ) ); 53 init( fileName.c_str( ) ); 54 } 54 } 55 /* 55 /* 56 ********************************************** 56 *************************************************************** 57 */ 57 */ 58 void G4GIDI_target::init( char const *fileName 58 void G4GIDI_target::init( char const *fileName ) { 59 59 60 int i, j, n, *p, *q, ir; 60 int i, j, n, *p, *q, ir; 61 MCGIDI_reaction *reaction; 61 MCGIDI_reaction *reaction; 62 62 63 smr_initialize( &smr, smr_status_Ok, 1 ); 63 smr_initialize( &smr, smr_status_Ok, 1 ); 64 sourceFilename = fileName; 64 sourceFilename = fileName; 65 target = MCGIDI_target_newRead( &smr, file 65 target = MCGIDI_target_newRead( &smr, fileName ); 66 if( !smr_isOk( &smr ) ) { 66 if( !smr_isOk( &smr ) ) { 67 smr_print( &smr, 1 ); 67 smr_print( &smr, 1 ); 68 throw 1; 68 throw 1; 69 } 69 } 70 projectilesPOPID = target->projectilePOP-> 70 projectilesPOPID = target->projectilePOP->globalPoPsIndex; 71 name = target->targetPOP->name; 71 name = target->targetPOP->name; 72 mass = G4GIDI_targetMass( target->targetPO 72 mass = G4GIDI_targetMass( target->targetPOP->name ); 73 equalProbableBinSampleMethod = "constant"; 73 equalProbableBinSampleMethod = "constant"; 74 elasticIndices = NULL; 74 elasticIndices = NULL; 75 nElasticIndices = nCaptureIndices = nFissi 75 nElasticIndices = nCaptureIndices = nFissionIndices = nOthersIndices = 0; 76 76 77 if( ( n = MCGIDI_target_numberOfReactions( 77 if( ( n = MCGIDI_target_numberOfReactions( &smr, target ) ) > 0 ) { 78 if( ( p = elasticIndices = (int *) smr 78 if( ( p = elasticIndices = (int *) smr_malloc2( &smr, n * sizeof( double ), 1, "elasticIndices" ) ) == NULL ) { 79 smr_print( &smr, 1 ); 79 smr_print( &smr, 1 ); 80 throw 1; 80 throw 1; 81 } 81 } 82 for( i = 0; i < n; i++ ) { /* Fin 82 for( i = 0; i < n; i++ ) { /* Find elastic channel(s). */ 83 reaction = MCGIDI_target_heated_ge 83 reaction = MCGIDI_target_heated_getReactionAtIndex( target->baseHeatedTarget, i ); 84 if( MCGIDI_reaction_getENDF_MTNumb 84 if( MCGIDI_reaction_getENDF_MTNumber( reaction ) == 2 ) { 85 *(p++) = i; 85 *(p++) = i; 86 nElasticIndices++; 86 nElasticIndices++; 87 } 87 } 88 } 88 } 89 captureIndices = p; 89 captureIndices = p; 90 for( i = 0; i < n; i++ ) { /* Fin 90 for( i = 0; i < n; i++ ) { /* Find capture channel(s). */ 91 reaction = MCGIDI_target_heated_ge 91 reaction = MCGIDI_target_heated_getReactionAtIndex( target->baseHeatedTarget, i ); 92 if( MCGIDI_reaction_getENDF_MTNumb 92 if( MCGIDI_reaction_getENDF_MTNumber( reaction ) == 102 ) { 93 *(p++) = i; 93 *(p++) = i; 94 nCaptureIndices++; 94 nCaptureIndices++; 95 } 95 } 96 } 96 } 97 97 98 fissionIndices = p; 98 fissionIndices = p; 99 for( i = 0; i < n; i++ ) { /* Fin 99 for( i = 0; i < n; i++ ) { /* Find fission channel(s). */ 100 reaction = MCGIDI_target_heated_ge 100 reaction = MCGIDI_target_heated_getReactionAtIndex( target->baseHeatedTarget, i ); 101 ir = MCGIDI_reaction_getENDF_MTNum 101 ir = MCGIDI_reaction_getENDF_MTNumber( reaction ); 102 if( ( ir != 18 ) && ( ir != 19 ) & 102 if( ( ir != 18 ) && ( ir != 19 ) && ( ir != 20 ) && ( ir != 21 ) && ( ir != 38 ) ) continue; 103 *(p++) = i; 103 *(p++) = i; 104 nFissionIndices++; 104 nFissionIndices++; 105 } 105 } 106 othersIndices = p; 106 othersIndices = p; 107 for( i = 0; i < n; i++ ) { /* Fin 107 for( i = 0; i < n; i++ ) { /* Find other channel(s). */ 108 for( j = 0, q = elasticIndices; j 108 for( j = 0, q = elasticIndices; j < nElasticIndices; j++, q++ ) if( *q == i ) break; 109 if( j < nElasticIndices ) continue 109 if( j < nElasticIndices ) continue; 110 for( j = 0, q = captureIndices; j 110 for( j = 0, q = captureIndices; j < nCaptureIndices; j++, q++ ) if( *q == i ) break; 111 if( j < nCaptureIndices ) continue 111 if( j < nCaptureIndices ) continue; 112 for( j = 0, q = fissionIndices; j 112 for( j = 0, q = fissionIndices; j < nFissionIndices; j++, q++ ) if( *q == i ) break; 113 if( j < nFissionIndices ) continue 113 if( j < nFissionIndices ) continue; 114 *p = i; 114 *p = i; 115 p++; 115 p++; 116 nOthersIndices++; 116 nOthersIndices++; 117 } 117 } 118 #if 0 118 #if 0 119 printf( "elastic %d: ", nElasticIndices ); 119 printf( "elastic %d: ", nElasticIndices ); 120 for( i = 0; i < nElasticIndices; i++ ) printf( 120 for( i = 0; i < nElasticIndices; i++ ) printf( " %d", elasticIndices[i] ); 121 printf( "\ncapture %d: ", nCaptureIndices ); 121 printf( "\ncapture %d: ", nCaptureIndices ); 122 for( i = 0; i < nCaptureIndices; i++ ) printf( 122 for( i = 0; i < nCaptureIndices; i++ ) printf( " %d", captureIndices[i] ); 123 printf( "\nfission %d: ", nFissionIndices ); 123 printf( "\nfission %d: ", nFissionIndices ); 124 for( i = 0; i < nFissionIndices; i++ ) printf( 124 for( i = 0; i < nFissionIndices; i++ ) printf( " %d", fissionIndices[i] ); 125 printf( "\nothers %d: ", nOthersIndices ); 125 printf( "\nothers %d: ", nOthersIndices ); 126 for( i = 0; i < nOthersIndices; i++ ) printf( 126 for( i = 0; i < nOthersIndices; i++ ) printf( " %d", othersIndices[i] ); 127 printf( "\n" ); 127 printf( "\n" ); 128 #endif 128 #endif 129 } 129 } 130 } 130 } 131 /* 131 /* 132 ********************************************** 132 *************************************************************** 133 */ 133 */ 134 G4GIDI_target::~G4GIDI_target( ) { 134 G4GIDI_target::~G4GIDI_target( ) { 135 135 136 MCGIDI_target_free( &smr, target ); 136 MCGIDI_target_free( &smr, target ); 137 smr_freeMemory( (void **) &elasticIndices 137 smr_freeMemory( (void **) &elasticIndices ); 138 smr_release( &smr ); 138 smr_release( &smr ); 139 } 139 } 140 /* 140 /* 141 ********************************************** 141 *************************************************************** 142 */ 142 */ 143 string *G4GIDI_target::getName( void ) { retur 143 string *G4GIDI_target::getName( void ) { return( &name ); } 144 /* 144 /* 145 ********************************************** 145 *************************************************************** 146 */ 146 */ 147 string *G4GIDI_target::getFilename( void ) { r 147 string *G4GIDI_target::getFilename( void ) { return( &sourceFilename ); } 148 /* 148 /* 149 ********************************************** 149 *************************************************************** 150 */ 150 */ 151 int G4GIDI_target::getZ( void ) { 151 int G4GIDI_target::getZ( void ) { 152 152 153 return( target->targetPOP->Z ); 153 return( target->targetPOP->Z ); 154 } 154 } 155 /* 155 /* 156 ********************************************** 156 *************************************************************** 157 */ 157 */ 158 int G4GIDI_target::getA( void ) { 158 int G4GIDI_target::getA( void ) { 159 159 160 return( target->targetPOP->A ); 160 return( target->targetPOP->A ); 161 } 161 } 162 /* 162 /* 163 ********************************************** 163 *************************************************************** 164 */ 164 */ 165 int G4GIDI_target::getM( void ) { 165 int G4GIDI_target::getM( void ) { 166 166 167 return( target->targetPOP->m ); 167 return( target->targetPOP->m ); 168 } 168 } 169 /* 169 /* 170 ********************************************** 170 *************************************************************** 171 */ 171 */ 172 double G4GIDI_target::getMass( void ) { 172 double G4GIDI_target::getMass( void ) { 173 173 174 return( mass ); 174 return( mass ); 175 } 175 } 176 /* 176 /* 177 ********************************************** 177 *************************************************************** 178 */ 178 */ 179 int G4GIDI_target::getTemperatures( double *te 179 int G4GIDI_target::getTemperatures( double *temperatures ) { 180 180 181 return( MCGIDI_target_getTemperatures( &sm 181 return( MCGIDI_target_getTemperatures( &smr, target, temperatures ) ); 182 } 182 } 183 /* 183 /* 184 ********************************************** 184 *************************************************************** 185 */ 185 */ 186 int G4GIDI_target::readTemperature( int index 186 int G4GIDI_target::readTemperature( int index ) { 187 187 188 return( MCGIDI_target_readHeatedTarget( &s 188 return( MCGIDI_target_readHeatedTarget( &smr, target, index ) ); 189 } 189 } 190 /* 190 /* 191 ********************************************** 191 *************************************************************** 192 */ 192 */ 193 string G4GIDI_target::getEqualProbableBinSampl 193 string G4GIDI_target::getEqualProbableBinSampleMethod( void ) { 194 194 195 return( equalProbableBinSampleMethod ); 195 return( equalProbableBinSampleMethod ); 196 } 196 } 197 /* 197 /* 198 ********************************************** 198 *************************************************************** 199 */ 199 */ 200 int G4GIDI_target::setEqualProbableBinSampleMe 200 int G4GIDI_target::setEqualProbableBinSampleMethod( string method ) { 201 201 202 if( method == "constant" ) { 202 if( method == "constant" ) { 203 equalProbableBinSampleMethod = "consta 203 equalProbableBinSampleMethod = "constant"; } 204 if( method == "linear" ) { 204 if( method == "linear" ) { 205 equalProbableBinSampleMethod = "linear 205 equalProbableBinSampleMethod = "linear"; } 206 else { 206 else { 207 return( 1 ); 207 return( 1 ); 208 } 208 } 209 return( 0 ); 209 return( 0 ); 210 } 210 } 211 /* 211 /* 212 ********************************************** 212 *************************************************************** 213 */ 213 */ 214 int G4GIDI_target::getNumberOfChannels( void ) 214 int G4GIDI_target::getNumberOfChannels( void ) { 215 215 216 return( MCGIDI_target_numberOfReactions( & 216 return( MCGIDI_target_numberOfReactions( &smr, target ) ); 217 } 217 } 218 /* 218 /* 219 ********************************************** 219 *************************************************************** 220 */ 220 */ 221 int G4GIDI_target::getNumberOfProductionChanne 221 int G4GIDI_target::getNumberOfProductionChannels( void ) { 222 222 223 return( MCGIDI_target_numberOfProductionRe 223 return( MCGIDI_target_numberOfProductionReactions( &smr, target ) ); 224 } 224 } 225 /* 225 /* 226 ********************************************** 226 *************************************************************** 227 */ 227 */ 228 channelID G4GIDI_target::getChannelsID( int ch 228 channelID G4GIDI_target::getChannelsID( int channelIndex ) { 229 229 230 MCGIDI_reaction *reaction; 230 MCGIDI_reaction *reaction; 231 231 232 if( ( reaction = MCGIDI_target_heated_getR 232 if( ( reaction = MCGIDI_target_heated_getReactionAtIndex_smr( &smr, target->baseHeatedTarget, channelIndex ) ) == NULL ) { 233 smr_print( &smr, 1 ); 233 smr_print( &smr, 1 ); 234 throw 1; 234 throw 1; 235 } 235 } 236 return( string( reaction->outputChannelStr 236 return( string( reaction->outputChannelStr ) ); /* Only works because channelID is defined to be string. */ 237 } 237 } 238 /* 238 /* 239 ********************************************** 239 *************************************************************** 240 */ 240 */ 241 vector<channelID> *G4GIDI_target::getChannelID 241 vector<channelID> *G4GIDI_target::getChannelIDs( void ) { 242 242 243 int i, n = MCGIDI_target_numberOfReactions 243 int i, n = MCGIDI_target_numberOfReactions( &smr, target ); 244 MCGIDI_reaction *reaction; 244 MCGIDI_reaction *reaction; 245 vector<channelID> *listOfChannels; 245 vector<channelID> *listOfChannels; 246 246 247 listOfChannels = new vector<channelID>( n 247 listOfChannels = new vector<channelID>( n ); 248 for( i = 0; i < n; i++ ) { 248 for( i = 0; i < n; i++ ) { 249 reaction = MCGIDI_target_heated_getRea 249 reaction = MCGIDI_target_heated_getReactionAtIndex( target->baseHeatedTarget, i ); 250 (*listOfChannels)[i] = reaction->outpu 250 (*listOfChannels)[i] = reaction->outputChannelStr; 251 } 251 } 252 return( listOfChannels ); 252 return( listOfChannels ); 253 } 253 } 254 /* 254 /* 255 ********************************************** 255 *************************************************************** 256 */ 256 */ 257 vector<channelID> *G4GIDI_target::getProductio 257 vector<channelID> *G4GIDI_target::getProductionChannelIDs( void ) { 258 258 259 return( NULL ); 259 return( NULL ); 260 } 260 } 261 /* 261 /* 262 ********************************************** 262 *************************************************************** 263 */ 263 */ 264 double G4GIDI_target::getTotalCrossSectionAtE( 264 double G4GIDI_target::getTotalCrossSectionAtE( double e_in, double temperature ) { 265 265 266 MCGIDI_quantitiesLookupModes mode( project 266 MCGIDI_quantitiesLookupModes mode( projectilesPOPID ); 267 267 268 mode.setProjectileEnergy( e_in ); 268 mode.setProjectileEnergy( e_in ); 269 mode.setCrossSectionMode( MCGIDI_quantityL 269 mode.setCrossSectionMode( MCGIDI_quantityLookupMode_pointwise ); 270 mode.setTemperature( temperature ); 270 mode.setTemperature( temperature ); 271 271 272 return( MCGIDI_target_getTotalCrossSection 272 return( MCGIDI_target_getTotalCrossSectionAtTAndE( NULL, target, mode, true ) ); 273 } 273 } 274 /* 274 /* 275 ********************************************** 275 *************************************************************** 276 */ 276 */ 277 double G4GIDI_target::getElasticCrossSectionAt 277 double G4GIDI_target::getElasticCrossSectionAtE( double e_in, double temperature ) { 278 278 279 return( sumChannelCrossSectionAtE( nElasti 279 return( sumChannelCrossSectionAtE( nElasticIndices, elasticIndices, e_in, temperature ) ); 280 } 280 } 281 /* 281 /* 282 ********************************************** 282 *************************************************************** 283 */ 283 */ 284 double G4GIDI_target::getCaptureCrossSectionAt 284 double G4GIDI_target::getCaptureCrossSectionAtE( double e_in, double temperature ) { 285 285 286 return( sumChannelCrossSectionAtE( nCaptur 286 return( sumChannelCrossSectionAtE( nCaptureIndices, captureIndices, e_in, temperature ) ); 287 } 287 } 288 /* 288 /* 289 ********************************************** 289 *************************************************************** 290 */ 290 */ 291 double G4GIDI_target::getFissionCrossSectionAt 291 double G4GIDI_target::getFissionCrossSectionAtE( double e_in, double temperature ) { 292 292 293 return( sumChannelCrossSectionAtE( nFissio 293 return( sumChannelCrossSectionAtE( nFissionIndices, fissionIndices, e_in, temperature ) ); 294 } 294 } 295 /* 295 /* 296 ********************************************** 296 *************************************************************** 297 */ 297 */ 298 double G4GIDI_target::getOthersCrossSectionAtE 298 double G4GIDI_target::getOthersCrossSectionAtE( double e_in, double temperature ) { 299 299 300 return( sumChannelCrossSectionAtE( nOthers 300 return( sumChannelCrossSectionAtE( nOthersIndices, othersIndices, e_in, temperature ) ); 301 } 301 } 302 /* 302 /* 303 ********************************************** 303 *************************************************************** 304 */ 304 */ 305 double G4GIDI_target::sumChannelCrossSectionAt 305 double G4GIDI_target::sumChannelCrossSectionAtE( int nIndices, int *indices, double e_in, double temperature ) { 306 306 307 int i; 307 int i; 308 double xsec = 0.; 308 double xsec = 0.; 309 MCGIDI_quantitiesLookupModes mode( project 309 MCGIDI_quantitiesLookupModes mode( projectilesPOPID ); 310 310 311 mode.setProjectileEnergy( e_in ); 311 mode.setProjectileEnergy( e_in ); 312 mode.setCrossSectionMode( MCGIDI_quantityL 312 mode.setCrossSectionMode( MCGIDI_quantityLookupMode_pointwise ); 313 mode.setTemperature( temperature ); 313 mode.setTemperature( temperature ); 314 314 315 for( i = 0; i < nIndices; i++ ) 315 for( i = 0; i < nIndices; i++ ) 316 xsec += MCGIDI_target_getIndexReaction 316 xsec += MCGIDI_target_getIndexReactionCrossSectionAtE( &smr, target, indices[i], mode, true ); 317 return( xsec ); 317 return( xsec ); 318 } 318 } 319 /* 319 /* 320 ********************************************** 320 *************************************************************** 321 */ 321 */ 322 int G4GIDI_target::sampleChannelCrossSectionAt 322 int G4GIDI_target::sampleChannelCrossSectionAtE( int nIndices, int *indices, double e_in, double temperature, 323 double (*rng)( void * ), void *rngStat 323 double (*rng)( void * ), void *rngState ) { 324 324 325 int i; 325 int i; 326 double xsec = 0., rxsec = sumChannelCrossS 326 double xsec = 0., rxsec = sumChannelCrossSectionAtE( nIndices, indices, e_in, temperature ) * rng( rngState ); 327 MCGIDI_quantitiesLookupModes mode( project 327 MCGIDI_quantitiesLookupModes mode( projectilesPOPID ); 328 328 329 mode.setProjectileEnergy( e_in ); 329 mode.setProjectileEnergy( e_in ); 330 mode.setCrossSectionMode( MCGIDI_quantityL 330 mode.setCrossSectionMode( MCGIDI_quantityLookupMode_pointwise ); 331 mode.setTemperature( temperature ); 331 mode.setTemperature( temperature ); 332 332 333 for( i = 0; i < nIndices - 1; i++ ) { 333 for( i = 0; i < nIndices - 1; i++ ) { 334 xsec += MCGIDI_target_getIndexReaction 334 xsec += MCGIDI_target_getIndexReactionCrossSectionAtE( &smr, target, indices[i], mode, true ); 335 if( xsec >= rxsec ) break; 335 if( xsec >= rxsec ) break; 336 } 336 } 337 return( indices[i] ); 337 return( indices[i] ); 338 } 338 } 339 /* 339 /* 340 ********************************************** 340 *************************************************************** 341 */ 341 */ 342 double G4GIDI_target::getElasticFinalState( do 342 double G4GIDI_target::getElasticFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) { 343 343 344 MCGIDI_decaySamplingInfo decaySamplingInfo 344 MCGIDI_decaySamplingInfo decaySamplingInfo; 345 MCGIDI_reaction *reaction = MCGIDI_target_ 345 MCGIDI_reaction *reaction = MCGIDI_target_heated_getReactionAtIndex_smr( &smr, target->baseHeatedTarget, elasticIndices[0] ); 346 MCGIDI_product *product; 346 MCGIDI_product *product; 347 MCGIDI_quantitiesLookupModes mode( project 347 MCGIDI_quantitiesLookupModes mode( projectilesPOPID ); 348 348 349 if( ( product = MCGIDI_outputChannel_getPr 349 if( ( product = MCGIDI_outputChannel_getProductAtIndex( &smr, &(reaction->outputChannel), 0 ) ) == NULL ) { 350 smr_print( &smr, 1 ); 350 smr_print( &smr, 1 ); 351 throw 1; 351 throw 1; 352 } 352 } 353 353 354 mode.setProjectileEnergy( e_in ); 354 mode.setProjectileEnergy( e_in ); 355 mode.setCrossSectionMode( MCGIDI_quantityL 355 mode.setCrossSectionMode( MCGIDI_quantityLookupMode_pointwise ); 356 mode.setTemperature( temperature ); 356 mode.setTemperature( temperature ); 357 357 358 decaySamplingInfo.isVelocity = 0; 358 decaySamplingInfo.isVelocity = 0; 359 decaySamplingInfo.rng = rng; 359 decaySamplingInfo.rng = rng; 360 decaySamplingInfo.rngState = rngState; 360 decaySamplingInfo.rngState = rngState; 361 if( MCGIDI_product_sampleMu( &smr, product 361 if( MCGIDI_product_sampleMu( &smr, product, mode, &decaySamplingInfo ) ) { 362 smr_print( &smr, 1 ); 362 smr_print( &smr, 1 ); 363 throw 1; 363 throw 1; 364 } 364 } 365 365 366 return( decaySamplingInfo.mu ); 366 return( decaySamplingInfo.mu ); 367 } 367 } 368 /* 368 /* 369 ********************************************** 369 *************************************************************** 370 */ 370 */ 371 vector<G4GIDI_Product> *G4GIDI_target::getCapt 371 vector<G4GIDI_Product> *G4GIDI_target::getCaptureFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) { 372 372 373 return( getFinalState( nCaptureIndices, ca 373 return( getFinalState( nCaptureIndices, captureIndices, e_in, temperature, rng, rngState ) ); 374 } 374 } 375 /* 375 /* 376 ********************************************** 376 *************************************************************** 377 */ 377 */ 378 vector<G4GIDI_Product> *G4GIDI_target::getFiss 378 vector<G4GIDI_Product> *G4GIDI_target::getFissionFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) { 379 379 380 return( getFinalState( nFissionIndices, fi 380 return( getFinalState( nFissionIndices, fissionIndices, e_in, temperature, rng, rngState ) ); 381 } 381 } 382 /* 382 /* 383 ********************************************** 383 *************************************************************** 384 */ 384 */ 385 vector<G4GIDI_Product> *G4GIDI_target::getOthe 385 vector<G4GIDI_Product> *G4GIDI_target::getOthersFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) { 386 386 387 return( getFinalState( nOthersIndices, oth 387 return( getFinalState( nOthersIndices, othersIndices, e_in, temperature, rng, rngState ) ); 388 } 388 } 389 /* 389 /* 390 ********************************************** 390 *************************************************************** 391 */ 391 */ 392 vector<G4GIDI_Product> *G4GIDI_target::getFina 392 vector<G4GIDI_Product> *G4GIDI_target::getFinalState( int nIndices, int *indices, double e_in, double temperature, 393 double (*rng)( void * ), void *rngState ) 393 double (*rng)( void * ), void *rngState ) { 394 394 395 int index = 0, i, n; 395 int index = 0, i, n; 396 vector<G4GIDI_Product> *products = NULL; 396 vector<G4GIDI_Product> *products = NULL; 397 MCGIDI_decaySamplingInfo decaySamplingInfo 397 MCGIDI_decaySamplingInfo decaySamplingInfo; 398 MCGIDI_sampledProductsDatas sampledProduct 398 MCGIDI_sampledProductsDatas sampledProductsDatas; 399 MCGIDI_sampledProductsData *productData; 399 MCGIDI_sampledProductsData *productData; 400 MCGIDI_quantitiesLookupModes mode( project 400 MCGIDI_quantitiesLookupModes mode( projectilesPOPID ); 401 401 402 decaySamplingInfo.isVelocity = 0; 402 decaySamplingInfo.isVelocity = 0; 403 decaySamplingInfo.rng = rng; 403 decaySamplingInfo.rng = rng; 404 decaySamplingInfo.rngState = rngState; 404 decaySamplingInfo.rngState = rngState; 405 405 406 if( nIndices == 0 ) { 406 if( nIndices == 0 ) { 407 return( NULL ); } 407 return( NULL ); } 408 else { 408 else { 409 if( nIndices == 1 ) { 409 if( nIndices == 1 ) { 410 index = indices[0]; } 410 index = indices[0]; } 411 else { 411 else { 412 index = sampleChannelCrossSectionA 412 index = sampleChannelCrossSectionAtE( nIndices, indices, e_in, temperature, rng, rngState ); 413 } 413 } 414 } 414 } 415 415 416 MCGIDI_sampledProducts_initialize( &smr, & 416 MCGIDI_sampledProducts_initialize( &smr, &sampledProductsDatas, 1000 ); 417 if( !smr_isOk( &smr ) ) { 417 if( !smr_isOk( &smr ) ) { 418 smr_print( &smr, 1 ); 418 smr_print( &smr, 1 ); 419 throw 1; 419 throw 1; 420 } 420 } 421 421 422 mode.setProjectileEnergy( e_in ); 422 mode.setProjectileEnergy( e_in ); 423 mode.setCrossSectionMode( MCGIDI_quantityL 423 mode.setCrossSectionMode( MCGIDI_quantityLookupMode_pointwise ); 424 mode.setTemperature( temperature ); 424 mode.setTemperature( temperature ); 425 425 426 n = MCGIDI_target_heated_sampleIndexReacti 426 n = MCGIDI_target_heated_sampleIndexReactionProductsAtE( &smr, target->baseHeatedTarget, index, mode, 427 &decaySamplingInfo, &sampledProduc 427 &decaySamplingInfo, &sampledProductsDatas ); 428 if( !smr_isOk( &smr ) ) { 428 if( !smr_isOk( &smr ) ) { 429 smr_print( &smr, 1 ); 429 smr_print( &smr, 1 ); 430 throw 1; 430 throw 1; 431 } 431 } 432 if( n > 0 ) { 432 if( n > 0 ) { 433 if( ( products = new vector<G4GIDI_Pro 433 if( ( products = new vector<G4GIDI_Product>( n ) ) != NULL ) { 434 for( i = 0; i < n; i++ ) { 434 for( i = 0; i < n; i++ ) { 435 productData = &(sampledProduct 435 productData = &(sampledProductsDatas.products[i]); 436 (*products)[i].A = productData 436 (*products)[i].A = productData->pop->A; 437 (*products)[i].Z = productData 437 (*products)[i].Z = productData->pop->Z; 438 (*products)[i].m = productData 438 (*products)[i].m = productData->pop->m; 439 (*products)[i].kineticEnergy = 439 (*products)[i].kineticEnergy = productData->kineticEnergy; 440 (*products)[i].px = productDat 440 (*products)[i].px = productData->px_vx; 441 (*products)[i].py = productDat 441 (*products)[i].py = productData->py_vy; 442 (*products)[i].pz = productDat 442 (*products)[i].pz = productData->pz_vz; 443 (*products)[i].birthTimeSec = produc 443 (*products)[i].birthTimeSec = productData->birthTimeSec; 444 } 444 } 445 } 445 } 446 } 446 } 447 MCGIDI_sampledProducts_release( &smr, &sam 447 MCGIDI_sampledProducts_release( &smr, &sampledProductsDatas ); 448 448 449 return( products ); 449 return( products ); 450 } 450 } 451 /* 451 /* 452 ********************************************** 452 *************************************************************** 453 */ 453 */ 454 double G4GIDI_target::getReactionsThreshold( i 454 double G4GIDI_target::getReactionsThreshold( int index ) { 455 455 456 return( MCGIDI_target_heated_getReactionsT 456 return( MCGIDI_target_heated_getReactionsThreshold( &smr, target->baseHeatedTarget, index ) ); 457 } 457 } 458 /* 458 /* 459 ********************************************** 459 *************************************************************** 460 */ 460 */ 461 double G4GIDI_target::getReactionsDomain( int 461 double G4GIDI_target::getReactionsDomain( int index, double *EMin, double *EMax ) { 462 462 463 return( MCGIDI_target_heated_getReactionsD 463 return( MCGIDI_target_heated_getReactionsDomain( &smr, target->baseHeatedTarget, index, EMin, EMax ) ); 464 } 464 } 465 465