Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lend/src/G4GIDI_target.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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