| 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