Geant4 Cross Reference |
1 #include <stdio.h> 1 2 #include <stdlib.h> 3 #include <string.h> 4 5 #include "PoPs.h" 6 #include "PoPs_private.h" 7 8 /* 9 In PoPs_addParticleIfNeeded and unitsDB_ad 10 lost if more memory cannot be allocated (n 11 */ 12 13 #if defined __cplusplus 14 namespace GIDI { 15 using namespace GIDI; 16 #endif 17 18 #define incrementalSize 1000 19 20 #define MeV2eV 1e6 21 #define MeV2keV 1e3 22 #define AMU2MeV 931.494028 23 #define AMU2eV ( MeV2eV * 931.494028 ) 24 #define K2MeV 8.6173856922566752e-11 25 #define K2eV ( MeV2eV * K2MeV ) 26 27 typedef struct unitConversions_s unitConversio 28 29 struct unitConversions_s { 30 char const *_from; 31 char const *_to; 32 double ratio; 33 }; 34 35 int PoPs_smr_ID = smr_unknownID; 36 static int referenceCount = 0; 37 static char versionStr[64] = ""; 38 39 /* 40 * For MPI the following need to be broadcast 41 */ 42 static unitsDB unitsRoot = { 0, 0, NULL }; 43 static PoPs popsRoot = { 0, 0, NULL, NULL }; 44 /* 45 * End need to MPI broadcasted. 46 */ 47 48 static unitConversions conversions[] = { { "am 49 { "MeV", "eV", MeV2eV }, { "MeV", "keV", M 50 51 static char const *PoPs_genreStrings[] = { "in 52 53 static int PoPs_particleProperIndex( int index 54 static int PoPs_sortedParticleIndex( char cons 55 static int unitsDB_release( void ); 56 /* 57 ============================================== 58 */ 59 const char *PoPs_version( void ) { 60 61 if( versionStr[0] == 0 ) snprintf( version 62 return( versionStr ); 63 } 64 /* 65 ============================================== 66 */ 67 int PoPs_versionMajor( void ) { return( POPS_V 68 int PoPs_versionMinor( void ) { return( POPS_V 69 int PoPs_versionPatchLevel( void ) { return( P 70 /* 71 ============================================== 72 */ 73 int PoPs_register( void ) { 74 75 if( referenceCount < 0 ) return( -1 ); 76 return( ++referenceCount ); 77 } 78 /* 79 ============================================== 80 */ 81 int PoPs_readDatabase( statusMessageReporting 82 83 return( PoPs_particleReadDatabase( smr, fi 84 } 85 /* 86 ============================================== 87 */ 88 int PoPs_release( statusMessageReporting *smr 89 90 referenceCount--; 91 if( referenceCount != 0 ) return( referenc 92 PoPs_releasePrivate( smr ); 93 return( 0 ); 94 } 95 /* 96 ============================================== 97 */ 98 int PoPs_releasePrivate( statusMessageReportin 99 100 int i; 101 102 for( i = 0; i < popsRoot.numberOfParticles 103 smr_freeMemory( (void **) &(popsRoot.pops) 104 popsRoot.sorted = NULL; 105 popsRoot.numberOfParticles = 0; 106 popsRoot.allocated = 0; 107 unitsDB_release( ); 108 return( 0 ); 109 } 110 /* 111 ============================================== 112 */ 113 PoP *PoPs_addParticleIfNeeded( statusMessageRe 114 /* 115 If particle with name pop->name is already 116 A NULL is returned if adding particle to p 117 */ 118 int i, index = PoPs_sortedParticleIndex( p 119 120 if( index >= 0 ) return( popsRoot.pops[PoP 121 if( popsRoot.numberOfParticles == popsRoot 122 int size = popsRoot.allocated + increm 123 PoP **sorted, **pops = (PoP **) smr_ma 124 125 if( pops == NULL ) return( NULL ); 126 sorted = &(pops[size]); 127 for( i = 0; i < popsRoot.numberOfParti 128 pops[i] = popsRoot.pops[i]; 129 sorted[i] = popsRoot.sorted[i]; 130 } 131 smr_freeMemory( (void **) &(popsRoot.p 132 popsRoot.pops = pops; 133 popsRoot.sorted = sorted; 134 popsRoot.allocated = size; 135 } 136 popsRoot.pops[popsRoot.numberOfParticles] 137 index = -index - 1; 138 for( i = popsRoot.numberOfParticles; i > i 139 popsRoot.sorted[index] = pop; 140 pop->index = popsRoot.numberOfParticles; 141 popsRoot.numberOfParticles++; 142 if( pop->genre == PoPs_genre_alias ) { /* 143 PoP *pop2; 144 145 for( pop2 = popsRoot.pops[pop->properI 146 pop2->aliasIndex = pop->index; 147 } 148 return( pop ); 149 } 150 /* 151 ============================================== 152 */ 153 PoP *PoPs_copyAddParticleIfNeeded( statusMessa 154 /* 155 If particle with name pop->name is already 156 If particle is not in popsRoot then copy p 157 A NULL is return if particle coping fails 158 */ 159 160 int index = PoPs_particleIndex( pop->name 161 PoP *newPoP; 162 163 if( index >= 0 ) return( popsRoot.pops[ind 164 165 if( ( newPoP = (PoP *) smr_malloc2( smr, s 166 if( PoP_copyParticle( smr, newPoP, pop ) ) 167 smr_freeMemory( (void **) &newPoP ); 168 return( NULL ); 169 } 170 if( PoPs_addParticleIfNeeded( smr, newPoP 171 PoP_free( newPoP ); 172 return( NULL ); 173 } 174 return( newPoP ); 175 } 176 /* 177 ============================================== 178 */ 179 PoP *PoPs_addAliasIfNeeded( statusMessageRepor 180 181 PoP *pop = PoP_makeAlias( smr, name, alias 182 183 if( pop != NULL ) { 184 if( pop->index < 0 ) { 185 if( PoPs_addParticleIfNeeded( smr, 186 PoP_free( pop ); 187 return( NULL ); 188 } 189 } 190 } 191 192 return( pop ); 193 } 194 /* 195 ============================================== 196 */ 197 int PoPs_numberOfParticle( void ) { 198 199 return( popsRoot.numberOfParticles ); 200 } 201 /* 202 ============================================== 203 */ 204 int PoPs_particleIndex( char const *name ) { 205 /* 206 A negative number is return if particle is 207 */ 208 int index = PoPs_sortedParticleIndex( name 209 210 if( index >= 0 ) index = PoPs_particleProp 211 return( index ); 212 } 213 /* 214 ============================================== 215 */ 216 int PoPs_particleIndex_smr( statusMessageRepor 217 218 int index = PoPs_particleIndex( name ); 219 220 if( index < 0 ) 221 smr_setReportError( smr, NULL, file, l 222 return( index ); 223 } 224 /* 225 ============================================== 226 */ 227 static int PoPs_particleProperIndex( int index 228 229 while( popsRoot.pops[index]->properIndex > 230 return( index ); 231 } 232 /* 233 ============================================== 234 */ 235 static int PoPs_sortedParticleIndex( char cons 236 /* 237 If name is a particle in popsRoot, its ind 238 a negative number is returned. For a parti 239 */ 240 int low = 0, mid, high = popsRoot.numberOf 241 242 if( high == 0 ) return( -1 ); 243 while( ( high - low ) > 1 ) { 244 mid = ( low + high ) >> 1; 245 iCmp = strcmp( name, popsRoot.sorted[m 246 if( iCmp == 0 ) return( mid ); 247 if( iCmp > 0 ) { 248 low = mid; } 249 else { 250 high = mid; 251 } 252 } // Loop checking, 11.05.2015, T. Koi 253 if( high == 1 ) { /* First point 254 if( !strcmp( name, popsRoot.sorted[0]- 255 if( strcmp( name, popsRoot.sorted[0]-> 256 } 257 if( high < popsRoot.numberOfParticles ) { 258 if( strcmp( name, popsRoot.sorted[high 259 } 260 return( -high - 1 ); 261 } 262 /* 263 ============================================== 264 */ 265 double PoPs_getMassInUnitOf( statusMessageRepo 266 267 int index = PoPs_particleIndex_smr( smr, n 268 269 if( index < 0 ) return( -1. ); 270 return( PoPs_getMassInUnitOf_atIndex( smr, 271 } 272 /* 273 ============================================== 274 */ 275 char const *PoPs_getName_atIndex( statusMessag 276 277 if( ( index < 0 ) || ( index >= popsRoot.n 278 smr_setReportError2( smr, PoPs_smr_ID, 279 return( NULL ); 280 } 281 return( popsRoot.pops[index]->name ); 282 } 283 /* 284 ============================================== 285 */ 286 double PoPs_getMassInUnitOf_atIndex( statusMes 287 288 double mass = -1.; 289 290 if( ( index < 0 ) || ( index >= popsRoot.n 291 smr_setReportError2( smr, PoPs_smr_ID, 292 else { 293 mass = PoP_getMassInUnitOf( smr, popsR 294 } 295 296 return( mass ); 297 } 298 /* 299 ============================================== 300 */ 301 enum PoPs_genre PoPs_getGenre( statusMessageRe 302 303 int index = PoPs_particleIndex_smr( smr, n 304 305 if( index < 0 ) return( PoPs_genre_invalid 306 return( popsRoot.pops[index]->genre ); 307 } 308 /* 309 ============================================== 310 */ 311 enum PoPs_genre PoPs_getGenre_atIndex( statusM 312 313 enum PoPs_genre genre = PoPs_genre_invalid 314 315 if( ( index < 0 ) || ( index >= popsRoot.n 316 smr_setReportError2( smr, PoPs_smr_ID, 317 else { 318 genre = popsRoot.pops[index]->genre; 319 } 320 return( genre ); 321 } 322 /* 323 ============================================== 324 */ 325 int PoPs_getZ_A_l( statusMessageReporting *smr 326 327 int index = PoPs_particleIndex_smr( smr, n 328 329 if( index < 0 ) return( -1 ); 330 return( PoPs_getZ_A_l_atIndex( smr, index, 331 } 332 /* 333 ============================================== 334 */ 335 int PoPs_getZ_A_l_atIndex( statusMessageReport 336 337 if( ( index < 0 ) || ( index >= popsRoot.n 338 smr_setReportError2( smr, PoPs_smr_ID, 339 return( -1 ); 340 } 341 *Z = popsRoot.pops[index]->Z; 342 *A = popsRoot.pops[index]->A; 343 *l = 0; 344 return( 0 ); 345 } 346 /* 347 ============================================== 348 */ 349 int PoPs_hasNucleus( statusMessageReporting *s 350 351 int index = PoPs_particleIndex_smr( smr, n 352 353 if( index < 0 ) return( -1 ); 354 return( PoPs_hasNucleus_atIndex( smr, inde 355 } 356 /* 357 ============================================== 358 */ 359 int PoPs_hasNucleus_atIndex( statusMessageRepo 360 /* 361 * If an error is encountered, a negative val 362 * contains a nucleus (is an atom, ion or nuc 363 */ 364 if( ( index < 0 ) || ( index >= popsRoot.n 365 smr_setReportError2( smr, PoPs_smr_ID, 366 return( -1 ); 367 } 368 if( ( popsRoot.pops[index]->genre == PoPs_ 369 if( protonIsNucleus ) { 370 if( strcmp( "p", popsRoot.pops[index]- 371 } 372 return( 0 ); 373 } 374 /* 375 ============================================== 376 */ 377 char const *PoPs_getAtomsName( statusMessageRe 378 379 int index = PoPs_particleIndex_smr( smr, n 380 381 if( index < 0 ) return( NULL ); 382 return( PoPs_getAtomsName_atIndex( smr, in 383 } 384 /* 385 ============================================== 386 */ 387 char const *PoPs_getAtomsName_atIndex( statusM 388 389 int atomIndex = PoPs_getAtomsIndex_atIndex 390 391 if( atomIndex < 0 ) return( NULL ); 392 return( popsRoot.pops[atomIndex]->name ); 393 } 394 /* 395 ============================================== 396 */ 397 int PoPs_getAtomsIndex( statusMessageReporting 398 399 int index = PoPs_particleIndex_smr( smr, n 400 401 if( index < 0 ) return( index ); 402 return( PoPs_getAtomsIndex_atIndex( smr, i 403 } 404 /* 405 ============================================== 406 */ 407 int PoPs_getAtomsIndex_atIndex( statusMessageR 408 409 char const *p = NULL; 410 411 if( ( index < 0 ) || ( index >= popsRoot.n 412 smr_setReportError2( smr, PoPs_smr_ID, 413 return( -1 ); 414 } 415 416 if( popsRoot.pops[index]->genre == PoPs_ge 417 418 if( strcmp( "p", popsRoot.pops[index] 419 p = "H1"; } 420 else { 421 if( popsRoot.pops[index]->genre != PoP 422 else if( strcmp( "h2", popsRoot.pops[i 423 p = "H2"; } 424 else if( strcmp( "h3", popsRoot.pops[i 425 p = "H3"; } 426 else if( strcmp( "he3", popsRoot.pops[ 427 p = "He3"; } 428 else if( strcmp( "he4", popsRoot.pops[ 429 p = "He4"; 430 } 431 } 432 if( p != NULL ) return( PoPs_particleIndex 433 return( -1 ); 434 } 435 /* 436 ============================================== 437 */ 438 PoP *PoPs_getParticle_atIndex( int index ) { 439 440 if( ( index < 0 ) || ( index >= popsRoot.n 441 return( popsRoot.pops[index] ); 442 } 443 /* 444 ============================================== 445 */ 446 char const *PoPs_genreTokenToString( enum PoPs 447 448 if( genre < PoPs_genre_invalid ) return( N 449 if( genre > PoPs_genre_atom ) return( NULL 450 return( PoPs_genreStrings[genre] ); 451 } 452 /* 453 ============================================== 454 */ 455 void PoPs_print( int sorted ) { 456 457 PoPs_write( stdout, sorted ); 458 } 459 /* 460 ============================================== 461 */ 462 void PoPs_write( FILE *f, int sorted ) { 463 464 int i1, properIndex; 465 PoP *pop; 466 467 fprintf( f, "Mass units: number of units = 468 for( i1 = 0; i1 < unitsRoot.numberOfUnits; 469 fprintf( f, " %s", unitsRoot.unsorted[ 470 } 471 fprintf( f, "\n\n" ); 472 fprintf( f, "Particles: number of particle 473 fprintf( f, " name in 474 fprintf( f, " 475 fprintf( f, " ---------------------------- 476 for( i1 = 0; i1 < popsRoot.numberOfParticl 477 if( sorted ) { 478 pop = popsRoot.sorted[i1]; } 479 else { 480 pop = popsRoot.pops[i1]; 481 } 482 properIndex = PoPs_particleProperIndex 483 fprintf( f, " %-24s %6d %-10s %15.8e 484 popsRoot.pops[properIndex]->mass, 485 if( PoPs_hasNucleus( NULL, pop->name, 486 fprintf( f, " T" ); } 487 else { 488 fprintf( f, " " ); 489 } 490 if( PoPs_hasNucleus( NULL, pop->name, 491 fprintf( f, " T" ); } 492 else { 493 fprintf( f, " " ); 494 } 495 if( pop->Z + pop->A > 0 ) { 496 fprintf( f, " %3d %3d", pop->Z, po 497 if( pop->l > 0 ) { 498 fprintf( f, " %d", pop->l ); } 499 else { 500 fprintf( f, " " ); 501 } } 502 else { 503 fprintf( f, " " ); 504 } 505 if( pop->genre == PoPs_genre_alias ) { 506 fprintf( f, " %s (%d)", popsRoot.p 507 else { 508 int aliasIndex; 509 510 for( aliasIndex = pop->aliasIndex; 511 } 512 fprintf( f, "\n" ); 513 } 514 } 515 516 /* 517 ========================== PoP functions 518 */ 519 /* 520 ============================================== 521 */ 522 PoP *PoP_new( statusMessageReporting *smr ) { 523 524 PoP *pop; 525 526 if( ( pop = (PoP *) smr_malloc2( smr, size 527 if( PoP_initialize( smr, pop ) != 0 ) pop 528 return( pop ); 529 } 530 /* 531 ============================================== 532 */ 533 int PoP_initialize( statusMessageReporting * / 534 535 pop->index = -1; 536 pop->properIndex = -1; 537 pop->aliasIndex = -1; 538 pop->genre = PoPs_genre_unknown; 539 pop->name = NULL; 540 pop->Z = 0; 541 pop->A = 0; 542 pop->mass = 0.0; 543 pop->massUnit = NULL; 544 return( 0 ); 545 } 546 /* 547 ============================================== 548 */ 549 int PoP_release( PoP *pop ) { 550 551 if( pop->name != NULL ) smr_freeMemory( (v 552 PoP_initialize( NULL, pop ); 553 return( 0 ); 554 } 555 /* 556 ============================================== 557 */ 558 PoP *PoP_free( PoP *pop ) { 559 560 PoP_release( pop ); 561 smr_freeMemory( (void **) &pop ); 562 return( NULL ); 563 } 564 /* 565 ============================================== 566 */ 567 int PoP_copyParticle( statusMessageReporting * 568 569 desc->index = -1; 570 desc->properIndex = src->properIndex; 571 desc->aliasIndex = src->aliasIndex; 572 desc->genre = src->genre; 573 if( ( desc->name = smr_allocateCopyString2 574 desc->Z = src->Z; 575 desc->A = src->A; 576 desc->l = src->l; 577 desc->mass = src->mass; 578 desc->massUnit = src->massUnit; 579 580 return( 0 ); 581 } 582 /* 583 ============================================== 584 */ 585 PoP *PoP_makeParticle( statusMessageReporting 586 587 PoP *pop; 588 589 if( ( pop = PoP_new( smr ) ) == NULL ) ret 590 if( ( pop->name = smr_allocateCopyString2( 591 PoP_free( pop ); 592 return( NULL ); 593 } 594 pop->genre = genre; 595 pop->mass = mass; 596 if( ( pop->massUnit = unitsDB_addUnitIfNee 597 return( pop ); 598 } 599 /* 600 ============================================== 601 */ 602 int PoP_setZ_A_l( statusMessageReporting * /*s 603 604 pop->Z = Z; 605 pop->A = A; 606 pop->l = l; 607 return( 0 ); 608 } 609 /* 610 ============================================== 611 */ 612 int PoP_getIndex( PoP *pop ) { 613 614 return( pop->index ); 615 } 616 /* 617 ============================================== 618 */ 619 char const *PoP_getName( PoP *pop ) { 620 621 return( pop->name ); 622 } 623 /* 624 ============================================== 625 */ 626 double PoP_getMassInUnitOf( statusMessageRepor 627 628 double mass = -1., ratio; 629 /* PoP *pop_ = pop;*/ 630 631 /*if( pop->genre == PoPs_genre_alias ) pop 632 if( PoPs_unitConversionRatio( pop->massUni 633 smr_setReportError2( smr, PoPs_smr_ID, 634 else { 635 mass = pop->mass * ratio; 636 } 637 638 return( mass ); 639 } 640 641 /* 642 ========================== alias functions 643 */ 644 /* 645 ============================================== 646 */ 647 PoP *PoP_makeAlias( statusMessageReporting *sm 648 649 int properIndex = PoPs_particleIndex( name 650 PoP *pop; 651 652 if( properIndex < 0 ) { 653 smr_setReportError2( smr, PoPs_smr_ID, 654 return( NULL ); 655 } 656 if( aliasIndex >= 0 ) { /* alias has a 657 PoP *truePop = popsRoot.pops[aliasInde 658 659 for( pop = truePop; strcmp( alias, pop 660 if( pop->genre != PoPs_genre_alias ) { 661 smr_setReportError2( smr, PoPs_smr 662 return( NULL ); 663 } 664 if( pop->properIndex != properIndex ) 665 smr_setReportError2( smr, PoPs_smr 666 alias, truePop->name, name ); 667 return( NULL ); 668 } } 669 else { 670 if( ( pop = PoP_new( smr ) ) == NULL ) 671 if( ( pop->name = smr_allocateCopyStri 672 PoP_free( pop ); 673 return( NULL ); 674 } 675 pop->properIndex = properIndex; 676 pop->genre = PoPs_genre_alias; 677 } 678 return( pop ); 679 } 680 681 /* 682 ========================== unitsDB functions 683 */ 684 /* 685 ============================================== 686 */ 687 static int unitsDB_release( void ) { 688 689 int i; 690 691 for( i = 0; i < unitsRoot.numberOfUnits; i 692 smr_freeMemory( (void **) &(unitsRoot.unso 693 unitsRoot.numberOfUnits = 0; 694 unitsRoot.allocated = 0; 695 return( 0 ); 696 } 697 /* 698 ============================================== 699 */ 700 char const *unitsDB_addUnitIfNeeded( statusMes 701 702 int i; 703 704 for( i = 0; i < unitsRoot.numberOfUnits; i 705 if( strcmp( unit, unitsRoot.unsorted[ 706 } 707 if( unitsRoot.numberOfUnits == unitsRoot.a 708 int size = unitsRoot.allocated + 20; 709 char const **unsorted = (char const ** 710 711 if( unsorted == NULL ) return( NULL ); 712 for( i = 0; i < unitsRoot.numberOfUnit 713 smr_freeMemory( (void **) &(unitsRoot. 714 unitsRoot.unsorted = unsorted; 715 unitsRoot.allocated = size; 716 } 717 if( ( unitsRoot.unsorted[unitsRoot.numberO 718 return( NULL ); 719 unitsRoot.numberOfUnits++; 720 return( unitsRoot.unsorted[unitsRoot.numbe 721 } 722 /* 723 ============================================== 724 */ 725 int unitsDB_index( statusMessageReporting * /* 726 727 int i; 728 729 for( i = 0; i < unitsRoot.numberOfUnits; i 730 if( !strcmp( unit, unitsRoot.unsorted[ 731 } 732 return( -1 ); 733 } 734 /* 735 ============================================== 736 */ 737 char const *unitsDB_stringFromIndex( statusMes 738 739 if( ( index < 0 ) || ( index >= unitsRoot. 740 smr_setReportError2( smr, PoPs_smr_ID, 741 return( NULL ); 742 } 743 return( unitsRoot.unsorted[index] ); 744 } 745 /* 746 ============================================== 747 */ 748 int PoPs_unitConversionRatio( char const *_fro 749 750 int i, n = sizeof( conversions ) / sizeof( 751 752 *ratio = 1.; 753 if( strcmp( _from, _to ) == 0 ) return( 0 754 for( i = 0; i < n; i++ ) { 755 if( strcmp( conversions[i]._from, _fro 756 if( strcmp( conversions[i]._to, _t 757 *ratio = conversions[i].ratio; 758 return( 0 ); 759 } } 760 else if( strcmp( conversions[i]._to, _ 761 if( strcmp( conversions[i]._from, 762 *ratio = 1. / conversions[i].r 763 return( 0 ); 764 } 765 } 766 } 767 return( 1 ); 768 } 769 #ifdef PoPs_MPI 770 #include "PoPs_Bcast_private.h" 771 /* 772 ============================================== 773 */ 774 int PoPs_Bcast( statusMessageReporting *smr, M 775 776 return( PoPs_Bcast2( smr, comm, rank, &uni 777 } 778 #endif 779 780 #if defined __cplusplus 781 } 782 #endif 783