00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062 #include <iostream>
00063 #include <stdlib.h>
00064
00065 #include "G4GIDI_target.hh"
00066 #include "G4GIDI_mass.hh"
00067 #include "G4GIDI_Misc.hh"
00068
00069 using namespace std;
00070 using namespace GIDI;
00071
00072
00073
00074
00075 G4GIDI_target::G4GIDI_target( const char *fileName ) {
00076
00077 init( fileName );
00078 }
00079
00080
00081
00082 G4GIDI_target::G4GIDI_target( string &fileName ) {
00083
00084 init( fileName.c_str( ) );
00085 }
00086
00087
00088
00089 void G4GIDI_target::init( const char *fileName ) {
00090
00091 int i, j, n, *p, *q;
00092 tpia_channel *channel;
00093
00094 smr_initialize( &smr );
00095 sourceFilename = fileName;
00096 target = tpia_target_createRead( &smr, fileName );
00097 if( !smr_isOk( &smr ) ) {
00098 smr_print( &smr, stderr, 1 );
00099 throw 1;
00100 }
00101 name = target->targetID->name;
00102 mass = G4GIDI_targetMass( target->targetID->name );
00103 equalProbableBinSampleMethod = "constant";
00104 elasticIndices = NULL;
00105 nElasticIndices = nCaptureIndices = nFissionIndices = nOthersIndices = 0;
00106
00107 if( ( n = tpia_target_numberOfChannels( &smr, target ) ) > 0 ) {
00108 p = elasticIndices = (int *) xData_malloc2( NULL, n * sizeof( double ), 1, "elasticIndices" );
00109 for( i = 0; i < n; i++ ) {
00110 channel = tpia_target_heated_getChannelAtIndex( target->baseHeatedTarget, i );
00111 if( channel->ENDL_C == 10 ) {
00112 *(p++) = i;
00113 nElasticIndices++;
00114 }
00115 }
00116 captureIndices = p;
00117 for( i = 0; i < n; i++ ) {
00118 channel = tpia_target_heated_getChannelAtIndex( target->baseHeatedTarget, i );
00119 if( channel->ENDL_C == 46 ) {
00120 *(p++) = i;
00121 nCaptureIndices++;
00122 }
00123 }
00124
00125 fissionIndices = p;
00126 for( i = 0; i < n; i++ ) {
00127 channel = tpia_target_heated_getChannelAtIndex( target->baseHeatedTarget, i );
00128 if( channel->fission != NULL ) {
00129 *(p++) = i;
00130 nFissionIndices++;
00131 }
00132 }
00133 othersIndices = p;
00134 for( i = 0; i < n; i++ ) {
00135 for( j = 0, q = elasticIndices; j < nElasticIndices; j++, q++ ) if( *q == i ) break;
00136 if( j < nElasticIndices ) continue;
00137 for( j = 0, q = captureIndices; j < nCaptureIndices; j++, q++ ) if( *q == i ) break;
00138 if( j < nCaptureIndices ) continue;
00139 for( j = 0, q = fissionIndices; j < nFissionIndices; j++, q++ ) if( *q == i ) break;
00140 if( j < nFissionIndices ) continue;
00141 *p = i;
00142 p++;
00143 nOthersIndices++;
00144 }
00145 }
00146 }
00147
00148
00149
00150 G4GIDI_target::~G4GIDI_target( ) {
00151
00152 tpia_target_free( &smr, target );
00153 xData_free( &smr, elasticIndices );
00154 smr_release( &smr );
00155 }
00156
00157
00158
00159 string *G4GIDI_target::getName( void ) { return( &name ); }
00160
00161
00162
00163 string *G4GIDI_target::getFilename( void ) { return( &sourceFilename ); }
00164
00165
00166
00167 int G4GIDI_target::getZ( void ) {
00168
00169 return( target->targetID->Z );
00170 }
00171
00172
00173
00174 int G4GIDI_target::getA( void ) {
00175
00176 return( target->targetID->A );
00177 }
00178
00179
00180
00181 int G4GIDI_target::getM( void ) {
00182
00183 return( target->targetID->m );
00184 }
00185
00186
00187
00188 double G4GIDI_target::getMass( void ) {
00189
00190 return( mass );
00191 }
00192
00193
00194
00195 int G4GIDI_target::getTemperatures( double *temperatures ) {
00196
00197 return( tpia_target_getTemperatures( &smr, target, temperatures ) );
00198 }
00199
00200
00201
00202 int G4GIDI_target::readTemperature( int index ) {
00203
00204 return( tpia_target_readHeatedTarget( &smr, target, index, 0 ) );
00205 }
00206
00207
00208
00209 string G4GIDI_target::getEqualProbableBinSampleMethod( void ) {
00210
00211 return( equalProbableBinSampleMethod );
00212 }
00213
00214
00215
00216 int G4GIDI_target::setEqualProbableBinSampleMethod( string method ) {
00217
00218 if( method == "constant" ) {
00219 equalProbableBinSampleMethod = "constant"; }
00220 if( method == "linear" ) {
00221 equalProbableBinSampleMethod = "linear"; }
00222 else {
00223 return( 1 );
00224 }
00225 return( 0 );
00226 }
00227
00228
00229
00230 int G4GIDI_target::getNumberOfChannels( void ) {
00231
00232 return( tpia_target_numberOfChannels( &smr, target ) );
00233 }
00234
00235
00236
00237 int G4GIDI_target::getNumberOfProductionChannels( void ) {
00238
00239 return( tpia_target_numberOfProductionChannels( &smr, target ) );
00240 }
00241
00242
00243
00244 vector<channelID> *G4GIDI_target::getChannelIDs( void ) {
00245
00246 return( getChannelIDs2( target->baseHeatedTarget->channels, tpia_target_numberOfChannels( &smr, target ) ) );
00247 }
00248
00249
00250
00251 vector<channelID> *G4GIDI_target::getProductionChannelIDs( void ) {
00252
00253 return( getChannelIDs2( target->baseHeatedTarget->productionChannels, tpia_target_numberOfProductionChannels( &smr, target ) ) );
00254 }
00255
00256
00257
00258 vector<channelID> *G4GIDI_target::getChannelIDs2( tpia_channel **channels, int n ) {
00259
00260 int i = 0;
00261 vector<channelID> *listOfChannels;
00262
00263 listOfChannels = new vector<channelID>( n );
00264 for( i = 0; i < n; i++ ) (*listOfChannels)[i].ID = channels[i]->outputChannel;
00265 return( listOfChannels );
00266 }
00267
00268
00269
00270 vector<double> *G4GIDI_target::getEnergyGridAtTIndex( int index ) {
00271
00272 xData_Int i, n;
00273 double *dEnergyGrid;
00274 vector<double> *energyGrid;
00275 vector<double>::iterator iter;
00276
00277 n = tpia_target_getEnergyGridAtTIndex( &smr, target, index, &dEnergyGrid );
00278 if( n < 0 ) return( NULL );
00279 energyGrid = new vector<double>( n );
00280 for( i = 0, iter = energyGrid->begin( ); i < n; i++, iter++ ) *iter = dEnergyGrid[i];
00281 return( energyGrid );
00282 }
00283
00284
00285
00286 double G4GIDI_target::getTotalCrossSectionAtE( double e_in, double temperature ) {
00287
00288 return( tpia_target_getTotalCrossSectionAtTAndE( NULL, target, temperature, -1, e_in, tpia_crossSectionType_pointwise ) );
00289 }
00290
00291
00292
00293 double G4GIDI_target::getElasticCrossSectionAtE( double e_in, double temperature ) {
00294
00295 return( sumChannelCrossSectionAtE( nElasticIndices, elasticIndices, e_in, temperature ) );
00296 }
00297
00298
00299
00300 double G4GIDI_target::getCaptureCrossSectionAtE( double e_in, double temperature ) {
00301
00302 return( sumChannelCrossSectionAtE( nCaptureIndices, captureIndices, e_in, temperature ) );
00303 }
00304
00305
00306
00307 double G4GIDI_target::getFissionCrossSectionAtE( double e_in, double temperature ) {
00308
00309 return( sumChannelCrossSectionAtE( nFissionIndices, fissionIndices, e_in, temperature ) );
00310 }
00311
00312
00313
00314 double G4GIDI_target::getOthersCrossSectionAtE( double e_in, double temperature ) {
00315
00316 return( sumChannelCrossSectionAtE( nOthersIndices, othersIndices, e_in, temperature ) );
00317 }
00318
00319
00320
00321 double G4GIDI_target::sumChannelCrossSectionAtE( int nIndices, int *indices, double e_in, double temperature ) {
00322
00323 int i;
00324 double xsec = 0.;
00325
00326 for( i = 0; i < nIndices; i++ )
00327 xsec += tpia_target_getIndexChannelCrossSectionAtE( &smr, target, indices[i], temperature, -1, e_in, tpia_crossSectionType_pointwise );
00328 return( xsec );
00329 }
00330
00331
00332
00333 int G4GIDI_target::sampleChannelCrossSectionAtE( int nIndices, int *indices, double e_in, double temperature,
00334 double (*rng)( void * ), void *rngState ) {
00335
00336 int i;
00337 double xsec = 0., rxsec = sumChannelCrossSectionAtE( nIndices, indices, e_in, temperature ) * tpia_misc_drng( rng, rngState );
00338
00339 for( i = 0; i < nIndices - 1; i++ ) {
00340 xsec += tpia_target_getIndexChannelCrossSectionAtE( &smr, target, indices[i], temperature, -1, e_in, tpia_crossSectionType_pointwise );
00341 if( xsec >= rxsec ) break;
00342 }
00343 return( indices[i] );
00344 }
00345
00346
00347
00348
00349 double G4GIDI_target::getElasticFinalState( double e_in, double , double (*rng)( void * ), void *rngState ) {
00350
00351 tpia_decaySamplingInfo decaySamplingInfo;
00352 tpia_channel *channel = tpia_target_heated_getChannelAtIndex_smr( &smr, target->baseHeatedTarget, elasticIndices[0] );
00353 tpia_product *product;
00354
00355 decaySamplingInfo.e_in = e_in;
00356 decaySamplingInfo.isVelocity = 0;
00357 tpia_frame_setColumn( &smr, &(decaySamplingInfo.frame), 0, tpia_referenceFrame_lab );
00358 decaySamplingInfo.samplingMethods = &(target->samplingMethods);
00359 decaySamplingInfo.rng = rng;
00360 decaySamplingInfo.rngState = rngState;
00361 product = tpia_decayChannel_getFirstProduct( &(channel->decayChannel) );
00362 tpia_angular_SampleMu( &smr, &(product->angular), &decaySamplingInfo );
00363
00364 return( decaySamplingInfo.mu );
00365 }
00366
00367
00368
00369 vector<G4GIDI_Product> *G4GIDI_target::getCaptureFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
00370
00371 return( getFinalState( nCaptureIndices, captureIndices, e_in, temperature, rng, rngState ) );
00372 }
00373
00374
00375
00376 vector<G4GIDI_Product> *G4GIDI_target::getFissionFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
00377
00378 return( getFinalState( nFissionIndices, fissionIndices, e_in, temperature, rng, rngState ) );
00379 }
00380
00381
00382
00383 vector<G4GIDI_Product> *G4GIDI_target::getOthersFinalState( double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
00384
00385 return( getFinalState( nOthersIndices, othersIndices, e_in, temperature, rng, rngState ) );
00386 }
00387
00388
00389
00390 vector<G4GIDI_Product> *G4GIDI_target::getFinalState( int nIndices, int *indices, double e_in, double temperature, double (*rng)( void * ), void *rngState ) {
00391
00392 #define nProductsMax 50
00393 int index = 0, i, n;
00394 vector<G4GIDI_Product> *products = NULL;
00395 tpia_decaySamplingInfo decaySamplingInfo;
00396 tpia_productOutgoingData productDatas[nProductsMax], *productData;
00397
00398 decaySamplingInfo.e_in = e_in;
00399 decaySamplingInfo.samplingMethods = &(target->samplingMethods);
00400 decaySamplingInfo.isVelocity = 0;
00401 tpia_frame_setColumn( &smr, &(decaySamplingInfo.frame), 0, tpia_referenceFrame_lab );
00402 decaySamplingInfo.rng = rng;
00403 decaySamplingInfo.rngState = rngState;
00404
00405 if( nIndices == 0 ) {
00406 return( NULL ); }
00407 else {
00408 if( nIndices == 1 ) {
00409 index = indices[0]; }
00410 else {
00411 index = sampleChannelCrossSectionAtE( nIndices, indices, e_in, temperature, rng, rngState );
00412 }
00413 }
00414 n = tpia_target_heated_sampleIndexChannelProductsAtE( &smr, target->baseHeatedTarget, index, &decaySamplingInfo, nProductsMax, productDatas );
00415 if( n > 0 ) {
00416 if( ( products = new vector<G4GIDI_Product>( n ) ) != NULL ) {
00417 for( i = 0; i < n; i++ ) {
00418 productData = &(productDatas[i]);
00419 (*products)[i].A = productData->productID->A;
00420 (*products)[i].Z = productData->productID->Z;
00421 (*products)[i].m = productData->productID->m;
00422 (*products)[i].kineticEnergy = productData->kineticEnergy;
00423 (*products)[i].px = productData->px_vx;
00424 (*products)[i].py = productData->py_vy;
00425 (*products)[i].pz = productData->pz_vz;
00426 }
00427 }
00428 }
00429
00430 return( products );
00431 #undef nProductsMax
00432 }