Geant4 Cross Reference |
1 /* 1 2 # <<BEGIN-copyright>> 3 # <<END-copyright>> 4 */ 5 #include <string.h> 6 #include <cmath> 7 8 #include "MCGIDI_fromTOM.h" 9 #include "MCGIDI_misc.h" 10 #include "MCGIDI_private.h" 11 12 #if defined __cplusplus 13 namespace GIDI { 14 using namespace GIDI; 15 #endif 16 17 /* 18 ********************************************** 19 */ 20 MCGIDI_angular *MCGIDI_angular_new( statusMess 21 22 MCGIDI_angular *angular; 23 24 if( ( angular = (MCGIDI_angular *) smr_mal 25 if( MCGIDI_angular_initialize( smr, angula 26 return( angular ); 27 } 28 /* 29 ********************************************** 30 */ 31 int MCGIDI_angular_initialize( statusMessageRe 32 33 memset( angular, 0, sizeof( MCGIDI_angular 34 return( 0 ); 35 } 36 /* 37 ********************************************** 38 */ 39 MCGIDI_angular *MCGIDI_angular_free( statusMes 40 41 MCGIDI_angular_release( smr, angular ); 42 smr_freeMemory( (void **) &angular ); 43 return( NULL ); 44 } 45 /* 46 ********************************************** 47 */ 48 int MCGIDI_angular_release( statusMessageRepor 49 50 51 MCGIDI_sampling_pdfsOfXGivenW_release( smr 52 53 MCGIDI_angular_initialize( smr, angular ); 54 return( 0 ); 55 } 56 /* 57 ********************************************** 58 */ 59 int MCGIDI_angular_setTwoBodyMasses( statusMes 60 double productMass_MeV, double residualMas 61 62 if( angular == NULL ) return( 0 ); 63 angular->projectileMass_MeV = projectileMa 64 angular->targetMass_MeV = targetMass_MeV; 65 angular->productMass_MeV = productMass_MeV 66 angular->residualMass_MeV = residualMass_M 67 return( 0 ); 68 } 69 /* 70 ********************************************** 71 */ 72 int MCGIDI_angular_parseFromTOM( statusMessage 73 74 MCGIDI_angular *angular = NULL; 75 xDataTOM_element *angularElement, *linearE 76 char const *nativeData; 77 ptwXYPoints *pdfXY = NULL; 78 ptwXPoints *cdfX = NULL; 79 ptwXY_interpolation interpolationXY, inter 80 81 if( ( angularElement = xDataTOME_getOneEle 82 if( ( angular = MCGIDI_angular_new( smr ) 83 84 if( ( nativeData = xDataTOM_getAttributesV 85 if( strcmp( nativeData, "isotropic" ) == 0 86 if( ( frameElement = xDataTOME_getOneE 87 smr_setReportError2( smr, smr_unkn 88 goto err; 89 } 90 angular->type = MCGIDI_angularType_iso 91 else if( strcmp( nativeData, "recoil" ) == 92 angular->type = MCGIDI_angularType_rec 93 else { 94 int i, j, n; 95 double norm, energyFactor; 96 nfu_status status; 97 xDataTOM_XYs *XYs; 98 xDataTOM_W_XYs *W_XYs; 99 ptwXYPoint *point; 100 MCGIDI_pdfsOfXGivenW *dists = &(angula 101 MCGIDI_pdfOfX *dist; 102 char const *energyUnit, *multiplicityP 103 104 if( ( linearElement = xDataTOME_getOne 105 if( ( linearElement = xDataTOME_ge 106 smr_setReportError2( smr, smr_ 107 goto err; 108 } 109 } 110 frameElement = linearElement; 111 112 if( MCGIDI_fromTOM_interpolation( smr, 113 if( MCGIDI_fromTOM_interpolation( smr, 114 dists->interpolationWY = interpolation 115 dists->interpolationXY = interpolation 116 117 if( ( W_XYs = (xDataTOM_W_XYs *) xData 118 if( ( dists->Ws = (double *) smr_mallo 119 if( ( dists->dist = (MCGIDI_pdfOfX *) 120 121 energyUnit = xDataTOM_subAxes_getUnit( 122 if( !smr_isOk( smr ) ) goto err; 123 energyFactor = MCGIDI_misc_getUnitConv 124 if( !smr_isOk( smr ) ) goto err; 125 126 for( i = 0; i < W_XYs->length; i++ ) { 127 XYs = &(W_XYs->XYs[i]); 128 dist = &(dists->dist[i]); 129 dists->Ws[i] = XYs->value * energy 130 if( ( pdfXY = MCGIDI_misc_dataFrom 131 if( ptwXY_simpleCoalescePoints( pd 132 dist->numberOfXs = n = (int) ptwXY 133 134 if( ( dist->Xs = (double *) smr_ma 135 dists->numberOfWs++; 136 dist->pdf = &(dist->Xs[n]); 137 dist->cdf = &(dist->pdf[n]); 138 139 for( j = 0; j < n; j++ ) { 140 point = ptwXY_getPointAtIndex_ 141 dist->Xs[j] = point->x; 142 dist->pdf[j] = point->y; 143 } 144 145 if( ( cdfX = ptwXY_runningIntegral 146 smr_setReportError2( smr, smr_ 147 goto err; 148 } 149 norm = ptwX_getPointAtIndex_Unsafe 150 if( norms != NULL ) { 151 ptwXY_setValueAtX( norms, XYs- 152 else if( std::fabs( 1. - norm ) > 153 smr_setReportError2( smr, smr_ 154 goto err; 155 } 156 for( j = 0; j < n; j++ ) dist->cdf 157 for( j = 0; j < n; j++ ) dist->pdf 158 pdfXY = ptwXY_free( pdfXY ); 159 cdfX = ptwX_free( cdfX ); 160 } 161 angular->type = MCGIDI_angularType_lin 162 } 163 164 if( frameElement != NULL ) { 165 if( ( angular->frame = MCGIDI_misc_get 166 } 167 168 distribution->angular = angular; 169 distribution->type = MCGIDI_distributionTy 170 171 return( 0 ); 172 173 err: 174 if( pdfXY != NULL ) ptwXY_free( pdfXY ); 175 if( cdfX != NULL ) cdfX = ptwX_free( cdfX 176 if ( angular != NULL ) MCGIDI_angular_free 177 return( 1 ); 178 } 179 /* 180 ********************************************** 181 */ 182 int MCGIDI_angular_sampleMu( statusMessageRepo 183 MCGIDI_decaySamplingInfo *decaySamplin 184 185 double randomMu = decaySamplingInfo->rng( 186 MCGIDI_pdfsOfXGivenW_sampled sampled; 187 188 switch( angular->type ) { 189 case MCGIDI_angularType_isotropic : 190 decaySamplingInfo->frame = angular->fr 191 decaySamplingInfo->mu = 1. - 2. * deca 192 break; 193 case MCGIDI_angularType_linear : 194 decaySamplingInfo->frame = angular->fr 195 sampled.smr = smr; 196 sampled.w = modes.getProjectileEnergy( 197 MCGIDI_sampling_sampleX_from_pdfsOfXGi 198 decaySamplingInfo->mu = sampled.x; 199 break; 200 case MCGIDI_angularType_recoil : 201 default : 202 smr_setReportError2( smr, smr_unknownI 203 } 204 return( !smr_isOk( smr ) ); 205 } 206 207 #if defined __cplusplus 208 } 209 #endif 210 211