ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MCGIDI_target.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MCGIDI_target.cc
1 /*
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5 
6 #include <map>
7 #include <string.h>
8 #include <cmath>
9 
10 #include "MCGIDI.h"
11 #include "MCGIDI_misc.h"
13 
14 #if defined __cplusplus
15  extern "C" {
16 namespace GIDI {
17 using namespace GIDI;
18 #endif
20 #if defined __cplusplus
21  }
22  }
23 #endif
24 /*
25 ************************************************************
26 */
27 
28 #if defined __cplusplus
29 namespace GIDI {
30 using namespace GIDI;
31 #endif
32 
34 
36 
37  if( ( target = (MCGIDI_target *) smr_malloc2( smr, sizeof( MCGIDI_target ), 0, "target" ) ) == NULL ) return( NULL );
38  if( MCGIDI_target_initialize( smr, target ) ) target = MCGIDI_target_free( smr, target );
39  return( target );
40 }
41 /*
42 ************************************************************
43 */
45 
46  memset( target, 0, sizeof( MCGIDI_target ) );
47  return( 0 );
48 }
49 /*
50 ************************************************************
51 */
53 
55 
56  if( ( target = MCGIDI_target_new( smr ) ) == NULL ) return( NULL );
57  if( MCGIDI_target_read( smr, target, fileName ) != 0 ) smr_freeMemory( (void **) &target );
58  return( target );
59 }
60 /*
61 ************************************************************
62 */
64  int projectile_PoPID, int target_PoPID ) {
65 
66  char *targetPath;
67 
68  if( ( targetPath = MCGIDI_map_findTargetViaPoPIDs( smr, map, evaluation, projectile_PoPID, target_PoPID ) ) == NULL ) return( 1 );
69  return( MCGIDI_target_read( smr, target, targetPath ) );
70 }
71 /*
72 ************************************************************
73 */
74 int MCGIDI_target_readFromMap( statusMessageReporting *smr, MCGIDI_target *target, MCGIDI_map *map, const char *evaluation, const char *projectileName,
75  const char *targetName ) {
76 
77  char *targetPath;
78 
79  if( ( targetPath = MCGIDI_map_findTarget( smr, map, evaluation, projectileName, targetName ) ) == NULL ) return( 1 );
80  return( MCGIDI_target_read( smr, target, targetPath ) );
81 }
82 /*
83 ************************************************************
84 */
86  int projectile_PoPID, int target_PoPID ) {
87 
88  char *targetPath;
90 
91  if( ( targetPath = MCGIDI_map_findTargetViaPoPIDs( smr, map, evaluation, projectile_PoPID, target_PoPID ) ) == NULL ) return( NULL );
92  target = MCGIDI_target_newRead( smr, targetPath );
93  smr_freeMemory( (void **) &targetPath );
94  return( target );
95 }
96 /*
97 ************************************************************
98 */
99 MCGIDI_target *MCGIDI_target_newReadFromMap( statusMessageReporting *smr, MCGIDI_map *map, const char *evaluation, const char *projectileName,
100  const char *targetName ) {
101 
102  char *targetPath;
104 
105  targetPath = MCGIDI_map_findTarget( smr, map, evaluation, projectileName, targetName );
106  if( targetPath == NULL ) return( NULL );
107  target = MCGIDI_target_newRead( smr, targetPath );
108  smr_freeMemory( (void **) &targetPath );
109  return( target );
110 }
111 /*
112 ************************************************************
113 */
115 
116  MCGIDI_target_release( smr, target );
117  smr_freeMemory( (void **) &target );
118  return( NULL );
119 }
120 /*
121 ************************************************************
122 */
124 
125  int i;
126 
127  smr_freeMemory( (void **) &(target->path) );
128  smr_freeMemory( (void **) &(target->absPath) );
129  xDataTOMAL_release( &(target->attributes) );
130  for( i = 0; i < target->nHeatedTargets; i++ ) {
131  smr_freeMemory( (void **) &(target->heatedTargets[i].path) );
132  smr_freeMemory( (void **) &(target->heatedTargets[i].contents) );
133  if( target->heatedTargets[i].heatedTarget != NULL ) MCGIDI_target_heated_free( smr, target->heatedTargets[i].heatedTarget );
134  }
135  smr_freeMemory( (void **) &(target->heatedTargets) );
136  smr_freeMemory( (void **) &(target->readHeatedTargets) );
137  MCGIDI_target_initialize( smr, target );
138  return( 0 );
139 }
140 /*
141 ************************************************************
142 */
143 int MCGIDI_target_read( statusMessageReporting *smr, MCGIDI_target *target, const char *fileName ) {
144 /*
145 * If a target has already been read into this target, user must have called MCGIDI_target_release before calling this routine.
146 * Otherwise, there will be memory leaks.
147 */
148  xDataXML_document *doc;
149  xDataXML_element *element, *child;
150  int i, iHeated, nHeated = 0, status = 1;
151  double temperature;
152  /* char *pReturnValue; */
153  char const *version, *contents;
154 
155  MCGIDI_target_initialize( smr, target );
156  if( ( target->path = smr_allocateCopyString2( smr, fileName, "path" ) ) == NULL ) return( status );
157  if( ( target->absPath = MCGIDI_misc_getAbsPath( smr, fileName ) ) == NULL ) return( _MCGIDI_target_releaseAndReturnOne( smr, target ) );
158  if( ( doc = xDataXML_importFile2( smr, fileName ) ) == NULL ) return( _MCGIDI_target_releaseAndReturnOne( smr, target ) );
159  element = xDataXML_getDocumentsElement( doc );
160  if( strcmp( element->name, "xTarget" ) != 0 ) {
161  MCGIDI_misc_setMessageError_Element( smr, NULL, element, __FILE__, __LINE__, 1, "input file's top element must be xTarget and not %s", element->name ); }
162  else {
163  status = 0;
164  if( ( version = xDataXML_getAttributesValueInElement( element, "version" ) ) == NULL ) {
165  smr_setReportError2( smr, smr_unknownID, 1, "version attribute missing from element '%s'", element->name );
166  status = 1; }
167  else {
168  if( strcmp( version, "xMCProcess 0.1" ) != 0 ) {
169  smr_setReportError2( smr, smr_unknownID, 1, "Unsupported version '%s' for element %s", version, element->name );
170  status = 1;
171  }
172  }
173  if( status == 0 ) {
174  /* pReturnValue = ( MCGIDI_misc_copyXMLAttributesToTOM( smr, &(target->attributes), &(element->attributes) ) ) ? NULL : target->path; */
175  for( nHeated = 0, child = xDataXML_getFirstElement( element ); child != NULL; nHeated++, child = xDataXML_getNextElement( child ) ) {
176  if( strcmp( child->name, "target" ) != 0 ) {
177  MCGIDI_misc_setMessageError_Element( smr, NULL, element, __FILE__, __LINE__, 1, "element can only have target sub-elements%s",
178  element->name );
179  status = 1;
180  break;
181  }
182  }
183  }
184  if( status == 0 ) {
185  if( ( target->heatedTargets = (MCGIDI_target_heated_info *) smr_malloc2( smr, nHeated * sizeof( MCGIDI_target_heated_info ), 1, "heatedTargets" ) ) == NULL ) {
186  status = 1; }
187  else {
188  if( ( target->readHeatedTargets = (MCGIDI_target_heated_info **) smr_malloc2( smr, nHeated * sizeof( MCGIDI_target_heated_info * ), 1, "heatedTargets" ) ) == NULL )
189  status = 1;
190  }
191  for( nHeated = 0, child = xDataXML_getFirstElement( element ); ( status == 0 ) && ( child != NULL );
192  nHeated++, child = xDataXML_getNextElement( child ) ) {
193  if( ( i = xDataXML_convertAttributeToDouble( smr, child, "temperature", &temperature, 1 ) ) != 0 ) {
194  if( i > 0 ) smr_setReportError2p( smr, smr_unknownID, 1, "target does not have a temperature attribute" );
195  status = 1;
196  break;
197  }
198  for( iHeated = 0; iHeated < nHeated; iHeated++ ) if( target->heatedTargets[iHeated].temperature > temperature ) break;
199  if( iHeated < nHeated ) for( i = nHeated; i >= iHeated; i-- ) target->heatedTargets[i+1] = target->heatedTargets[i];
200  target->heatedTargets[iHeated].temperature = temperature;
201  target->heatedTargets[iHeated].path = NULL;
202  target->heatedTargets[iHeated].contents = NULL;
203  target->heatedTargets[iHeated].heatedTarget = NULL;
204  if( ( contents = xDataXML_getAttributesValueInElement( child, "contents" ) ) != NULL ) {
205  if( ( target->heatedTargets[iHeated].contents = smr_allocateCopyString2( smr, contents, "contents" ) ) == NULL ) {
206  status = 1;
207  break;
208  }
209  }
210  if( ( contents = xDataXML_getAttributesValueInElement( child, "file" ) ) == NULL ) {
211  status = 1;
212  break;
213  }
214  if( ( target->heatedTargets[iHeated].path = (char *) smr_malloc2(smr, strlen( target->absPath ) + strlen( contents ) + 2, 0, "path") ) == NULL) {
215  status = 1;
216  break;
217  }
218  strcpy( target->heatedTargets[iHeated].path, target->absPath );
219  *strrchr( target->heatedTargets[iHeated].path, '/' ) = 0;
220  strcat( target->heatedTargets[iHeated].path, "/" );
221  strcat( target->heatedTargets[iHeated].path, contents );
222  target->nHeatedTargets++;
223  }
224  }
225  }
226  xDataXML_freeDoc( smr, doc );
227  if( status == 0 ) {
228  for( i = 0; i < nHeated; i++ ) target->heatedTargets[i].ordinal = i;
229  for( i = 0; i < nHeated; i++ ) if( target->heatedTargets[i].contents == NULL ) break;
230  if( i == nHeated ) i = 0; /* All heated targets are crossSection only. */
231  if( MCGIDI_target_readHeatedTarget( smr, target, i ) == 0 ) {
232  target->baseHeatedTarget = target->heatedTargets[i].heatedTarget; }
233  else {
234  MCGIDI_target_release( NULL, target );
235  status = 1;
236  } }
237  else {
238  MCGIDI_target_release( smr, target );
239  }
240  return( status );
241 }
242 /*
243 ************************************************************
244 */
246 
247  return( xDataTOMAL_getAttributesValue( &(target->attributes), name ) );
248 }
249 /*
250 ************************************************************
251 */
253 
254  int i;
255 
256  if( temperatures != NULL ) for( i = 0; i < target->nHeatedTargets; i++ ) temperatures[i] = target->heatedTargets[i].temperature;
257  return( target->nHeatedTargets );
258 }
259 /*
260 ************************************************************
261 */
263 
264  int i;
265 
266  if( ( index < 0 ) || ( index >= target->nHeatedTargets ) ) {
267  smr_setReportError2( smr, smr_unknownID, 1, "temperature index = %d out of range (0 <= index < %d", index, target->nHeatedTargets );
268  return( -1 );
269  }
270  if( target->heatedTargets[index].heatedTarget != NULL ) return( 1 );
271  if( ( target->heatedTargets[index].heatedTarget = MCGIDI_target_heated_newRead( smr, target->heatedTargets[index].path ) ) != NULL ) {
272  target->projectilePOP = target->heatedTargets[index].heatedTarget->projectilePOP;
273  target->targetPOP = target->heatedTargets[index].heatedTarget->targetPOP;
274  if( target->heatedTargets[index].heatedTarget != NULL ) {
275  target->heatedTargets[index].heatedTarget->ordinal = target->heatedTargets[index].ordinal;
276  for( i = target->nReadHeatedTargets; i > 0; i-- ) {
277  if( target->readHeatedTargets[i-1]->temperature < target->heatedTargets[index].temperature ) break;
278  target->readHeatedTargets[i] = target->readHeatedTargets[i-1];
279  }
280  target->readHeatedTargets[i] = &(target->heatedTargets[i]);
281  target->nReadHeatedTargets++;
282  }
283  }
284  return( ( target->heatedTargets[index].heatedTarget == NULL ? -1 : 0 ) );
285 }
286 /*
287 ************************************************************
288 */
290 
291  if( ( index < 0 ) || ( index >= target->nHeatedTargets ) ) {
292  smr_setReportError2( smr, smr_unknownID, 1, "temperature index = %d out of range (0 <= index < %d", index, target->nHeatedTargets );
293  return( NULL );
294  }
295  if( target->heatedTargets[index].heatedTarget == NULL ) MCGIDI_target_readHeatedTarget( smr, target, index );
296  return( target->heatedTargets[index].heatedTarget );
297 }
298 /*
299 ************************************************************
300 */
302 
303  if( ( index < 0 ) || ( index >= target->nHeatedTargets ) ) {
304  smr_setReportError2( smr, smr_unknownID, 1, "temperature index = %d out of range (0 <= index < %d", index, target->nHeatedTargets );
305  return( NULL );
306  }
307  if( target->heatedTargets[index].heatedTarget == NULL ) {
308  smr_setReportError2( smr, smr_unknownID, 1, "temperature index = %d not read in", index );
309  return( NULL );
310  }
311  return( target->heatedTargets[index].heatedTarget );
312 }
313 /*
314 ************************************************************
315 */
317 
319 }
320 /*
321 ************************************************************
322 */
324 
326 
327  if( reaction == NULL ) return( MCGIDI_reactionType_unknown_e );
328  return( MCGIDI_reaction_getReactionType( smr, reaction ) );
329 }
330 /*
331 ************************************************************
332 */
334 
335  return( MCGIDI_target_heated_getReactionAtIndex( target->baseHeatedTarget, index ) );
336 }
337 /*
338 ************************************************************
339 */
341 
342  return( MCGIDI_target_heated_getReactionAtIndex_smr( smr, target->baseHeatedTarget, index ) );
343 }
344 /*
345 ************************************************************
346 */
348 
349 #if 0
350  return( MCGIDI_target_heated_numberOfProductionReactions( smr, target->baseHeatedTarget ) );
351 #endif
352 return( 0 );
353 }
354 
355 /*
356 ************************************************************
357 */
359  bool sampling ) {
360 
361  int i;
362  double xsec = 0., xsec1, xsec2, temperature = modes.getTemperature( );
363 
364  for( i = 0; i < target->nReadHeatedTargets; i++ ) if( target->readHeatedTargets[i]->temperature > temperature ) break;
365  if( i == 0 ) {
366  xsec = MCGIDI_target_heated_getTotalCrossSectionAtE( smr, target->readHeatedTargets[0]->heatedTarget, modes, sampling ); }
367  else if( i == target->nReadHeatedTargets ) {
368  xsec = MCGIDI_target_heated_getTotalCrossSectionAtE( smr, target->readHeatedTargets[i-1]->heatedTarget, modes, sampling ); }
369  else {
370  xsec1 = MCGIDI_target_heated_getTotalCrossSectionAtE( smr, target->readHeatedTargets[i-1]->heatedTarget, modes, sampling );
371  xsec2 = MCGIDI_target_heated_getTotalCrossSectionAtE( smr, target->readHeatedTargets[i ]->heatedTarget, modes, sampling );
372  xsec = ( ( target->readHeatedTargets[i]->temperature - temperature ) * xsec1 +
373  ( temperature - target->readHeatedTargets[i-1]->temperature ) * xsec2 ) /
374  ( target->readHeatedTargets[i]->temperature - target->readHeatedTargets[i-1]->temperature );
375  }
376 
377  return( xsec );
378 }
379 /*
380 ************************************************************
381 */
382 int MCGIDI_target_getDomain( statusMessageReporting *smr, MCGIDI_target *target, double *EMin, double *EMax ) {
383 
384  int ir, nr = MCGIDI_target_numberOfReactions( smr, target );
385  double EMin_, EMax_;
386 
387  for( ir = 0; ir < nr; ir++ ) {
388  MCGIDI_target_heated_getReactionsDomain( smr, target->baseHeatedTarget, ir, &EMin_, &EMax_ );
389  if( ir == 0 ) {
390  *EMin = EMin_;
391  *EMax = EMax_; }
392  else {
393  if( *EMin > EMin_ ) *EMin = EMin_;
394  if( *EMax < EMax_ ) *EMax = EMax_;
395  }
396  }
397  return( 0 );
398 }
399 /*
400 ************************************************************
401 */
403  bool sampling ) {
404 
405  int i;
406  double xsec = 0., xsec1, xsec2, temperature = modes.getTemperature( );
407 
408  for( i = 0; i < target->nReadHeatedTargets; i++ ) if( target->readHeatedTargets[i]->temperature > temperature ) break;
409  if( i == 0 ) {
410  xsec = MCGIDI_target_heated_getIndexReactionCrossSectionAtE( smr, target->readHeatedTargets[0]->heatedTarget, index, modes, sampling ); }
411  else if( i == target->nReadHeatedTargets ) {
412  xsec = MCGIDI_target_heated_getIndexReactionCrossSectionAtE( smr, target->readHeatedTargets[i-1]->heatedTarget, index, modes, sampling ); }
413  else {
414  xsec1 = MCGIDI_target_heated_getIndexReactionCrossSectionAtE(smr, target->readHeatedTargets[i-1]->heatedTarget, index, modes, sampling );
415  xsec2 = MCGIDI_target_heated_getIndexReactionCrossSectionAtE(smr, target->readHeatedTargets[i ]->heatedTarget, index, modes, sampling );
416  xsec = ( ( target->readHeatedTargets[i]->temperature - temperature ) * xsec1 +
417  ( temperature - target->readHeatedTargets[i-1]->temperature ) * xsec2 ) /
418  ( target->readHeatedTargets[i]->temperature - target->readHeatedTargets[i-1]->temperature );
419  }
420 
421  return( xsec );
422 }
423 /*
424 ************************************************************
425 */
427  double (*userrng)( void * ), void *rngState ) {
428 
429  int ir, nr = MCGIDI_target_numberOfReactions( smr, target );
430  double rngValue = (*userrng)( rngState );
431  double cumm_xsec = 0., r_xsec = rngValue * totalXSec;
432 
433  for( ir = 0; ir < nr; ir++ ) {
434  cumm_xsec += MCGIDI_target_getIndexReactionCrossSectionAtE( smr, target, ir, modes, true );
435  if( cumm_xsec >= r_xsec ) break;
436  }
437  if( ir == nr ) {
438  if( ( totalXSec - cumm_xsec ) >= 1e-12 * totalXSec ) {
440  "Failed to sample a reaction for temperature = %.12e, energy = %.12e, totalXSec = %16.e, rngValue = %16.e, r_xsec = %16.e, cumm_xsec = %16.e",
441  modes.getTemperature( ), modes.getProjectileEnergy( ), totalXSec, rngValue, r_xsec, cumm_xsec );
442  return( -1 );
443  }
444  ir--; /* May not be correct but close. */
445  }
448 
449  if( modes.getGroupIndex( ) == reaction->thresholdGroupIndex ) {
450  double dEnergy = modes.getProjectileEnergy( ) - reaction->EMin;
451 
452  if( dEnergy <= 0 ) return( MCGIDI_nullReaction );
453  if( ( (*userrng)( rngState ) * reaction->thresholdGroupDomain ) > dEnergy ) return( MCGIDI_nullReaction );
454  }
455  }
456  return( ir );
457 }
458 /*
459 ************************************************************
460 */
462  MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo, MCGIDI_sampledProductsDatas *productDatas ) {
463 
464  MCGIDI_sampledProductsData productData;
465 
466  productData.isVelocity = decaySamplingInfo->isVelocity;
467  productData.pop = target->projectilePOP;
468  productData.kineticEnergy = modes.getProjectileEnergy( );
469  productData.px_vx = 0.;
470  productData.py_vy = 0.;
471  productData.pz_vz = std::sqrt( productData.kineticEnergy * ( productData.kineticEnergy + 2. * productData.pop->mass_MeV ) );
472  if( productData.isVelocity ) productData.pz_vz *=
473  MCGIDI_speedOfLight_cm_sec / std::sqrt( productData.pz_vz * productData.pz_vz + productData.pop->mass_MeV * productData.pop->mass_MeV );
474  productData.delayedNeutronIndex = 0;
475  productData.delayedNeutronRate = 0.;
476  productData.birthTimeSec = 0;
477 
478  productDatas->numberOfProducts = 0;
479  MCGIDI_sampledProducts_addProduct( smr, productDatas, &productData );
480  return( productDatas->numberOfProducts );
481 }
482 /*
483 ************************************************************
484 */
486  MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo, MCGIDI_sampledProductsDatas *productData ) {
487 
488  return( MCGIDI_target_heated_sampleIndexReactionProductsAtE( smr, target->baseHeatedTarget, index, modes, decaySamplingInfo, productData ) );
489 }
490 /*
491 ************************************************************
492 */
495 
496  return( MCGIDI_target_heated_getIndexReactionFinalQ( smr, target->baseHeatedTarget, index, modes ) );
497 }
498 /*
499 ************************************************************
500 */
501 std::map<int, enum MCGIDI_transportability> const *MCGIDI_target_getUniqueProducts( statusMessageReporting *smr, MCGIDI_target *target ) {
502 
504 }
505 /*
506 ************************************************************
507 */
509 
510  int i1, status = 0;
511 
512  for( i1 = 0; i1 < target->nReadHeatedTargets; i1++ ) {
513  if( ( status = MCGIDI_target_heated_recast( smr, target->readHeatedTargets[i1]->heatedTarget, settings ) ) != 0 ) break;
514  }
515  return( status );
516 }
517 /*
518 ************************************************************
519 */
521 
522  MCGIDI_target_release( smr, target );
523  return( 1 );
524 }
525 
526 #if defined __cplusplus
527 }
528 #endif