Geant4 Cross Reference |
1 /* 1 /* 2 # <<BEGIN-copyright>> 2 # <<BEGIN-copyright>> 3 # <<END-copyright>> 3 # <<END-copyright>> 4 */ 4 */ 5 #include <string.h> 5 #include <string.h> 6 6 7 #include "MCGIDI_fromTOM.h" 7 #include "MCGIDI_fromTOM.h" 8 #include "MCGIDI.h" 8 #include "MCGIDI.h" 9 #include "MCGIDI_misc.h" 9 #include "MCGIDI_misc.h" 10 10 11 #if defined __cplusplus 11 #if defined __cplusplus 12 namespace GIDI { 12 namespace GIDI { 13 using namespace GIDI; 13 using namespace GIDI; 14 #endif 14 #endif 15 15 16 static int MCGIDI_energyAngular_linear_parseFr 16 static int MCGIDI_energyAngular_linear_parseFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution ); 17 /* 17 /* 18 ********************************************** 18 ************************************************************ 19 */ 19 */ 20 MCGIDI_energyAngular *MCGIDI_energyAngular_new 20 MCGIDI_energyAngular *MCGIDI_energyAngular_new( statusMessageReporting *smr ) { 21 21 22 MCGIDI_energyAngular *energyAngular; 22 MCGIDI_energyAngular *energyAngular; 23 23 24 if( ( energyAngular = (MCGIDI_energyAngula 24 if( ( energyAngular = (MCGIDI_energyAngular *) smr_malloc2( smr, sizeof( MCGIDI_energyAngular ), 0, "energyAngular" ) ) == NULL ) return( NULL ); 25 if( MCGIDI_energyAngular_initialize( smr, 25 if( MCGIDI_energyAngular_initialize( smr, energyAngular ) ) energyAngular = MCGIDI_energyAngular_free( smr, energyAngular ); 26 return( energyAngular ); 26 return( energyAngular ); 27 } 27 } 28 /* 28 /* 29 ********************************************** 29 ************************************************************ 30 */ 30 */ 31 int MCGIDI_energyAngular_initialize( statusMes 31 int MCGIDI_energyAngular_initialize( statusMessageReporting * /*smr*/, MCGIDI_energyAngular *energyAngular ) { 32 32 33 memset( energyAngular, 0, sizeof( MCGIDI_e 33 memset( energyAngular, 0, sizeof( MCGIDI_energyAngular ) ); 34 return( 0 ); 34 return( 0 ); 35 } 35 } 36 /* 36 /* 37 ********************************************** 37 ************************************************************ 38 */ 38 */ 39 MCGIDI_energyAngular *MCGIDI_energyAngular_fre 39 MCGIDI_energyAngular *MCGIDI_energyAngular_free( statusMessageReporting *smr, MCGIDI_energyAngular *energyAngular ) { 40 40 41 MCGIDI_energyAngular_release( smr, energyA 41 MCGIDI_energyAngular_release( smr, energyAngular ); 42 smr_freeMemory( (void **) &energyAngular ) 42 smr_freeMemory( (void **) &energyAngular ); 43 return( NULL ); 43 return( NULL ); 44 } 44 } 45 /* 45 /* 46 ********************************************** 46 ************************************************************ 47 */ 47 */ 48 int MCGIDI_energyAngular_release( statusMessag 48 int MCGIDI_energyAngular_release( statusMessageReporting *smr, MCGIDI_energyAngular *energyAngular ) { 49 49 50 int i; 50 int i; 51 51 52 for( i = 0; i < energyAngular->pdfOfEpGive 52 for( i = 0; i < energyAngular->pdfOfEpGivenE.numberOfWs; i++ ) { 53 MCGIDI_sampling_pdfsOfXGivenW_release( 53 MCGIDI_sampling_pdfsOfXGivenW_release( smr, &(energyAngular->pdfOfMuGivenEAndEp[i]) ); 54 } 54 } 55 smr_freeMemory( (void **) &(energyAngular- 55 smr_freeMemory( (void **) &(energyAngular->pdfOfMuGivenEAndEp) ); 56 MCGIDI_sampling_pdfsOfXGivenW_release( smr 56 MCGIDI_sampling_pdfsOfXGivenW_release( smr, &(energyAngular->pdfOfEpGivenE) ); 57 MCGIDI_energyAngular_initialize( smr, ener 57 MCGIDI_energyAngular_initialize( smr, energyAngular ); 58 58 59 return( 0 ); 59 return( 0 ); 60 } 60 } 61 /* 61 /* 62 ********************************************** 62 ************************************************************ 63 */ 63 */ 64 int MCGIDI_energyAngular_parseFromTOM( statusM 64 int MCGIDI_energyAngular_parseFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution ) { 65 65 66 xDataTOM_element *energyAngularElement; 66 xDataTOM_element *energyAngularElement; 67 char const *nativeData; 67 char const *nativeData; 68 68 69 if( ( energyAngularElement = xDataTOME_get 69 if( ( energyAngularElement = xDataTOME_getOneElementByName( smr, element, "energyAngular", 1 ) ) == NULL ) goto err; 70 70 71 if( ( nativeData = xDataTOM_getAttributesV 71 if( ( nativeData = xDataTOM_getAttributesValueInElement( energyAngularElement, "nativeData" ) ) == NULL ) goto err; 72 if( strcmp( nativeData, "KalbachMann" ) == 72 if( strcmp( nativeData, "KalbachMann" ) == 0 ) { 73 return( MCGIDI_KalbachMann_parseFromTO 73 return( MCGIDI_KalbachMann_parseFromTOM( smr, energyAngularElement, distribution ) ); } 74 else if( strcmp( nativeData, "linear" ) == 74 else if( strcmp( nativeData, "linear" ) == 0 ) { 75 return( MCGIDI_energyAngular_linear_pa 75 return( MCGIDI_energyAngular_linear_parseFromTOM( smr, energyAngularElement, distribution ) ); } 76 else { 76 else { 77 smr_setReportError2( smr, smr_unknownI 77 smr_setReportError2( smr, smr_unknownID, 1, "energyAngular nativeData = '%s' not supported", nativeData ); 78 goto err; 78 goto err; 79 } 79 } 80 80 81 return( 0 ); 81 return( 0 ); 82 82 83 err: 83 err: 84 return( 1 ); 84 return( 1 ); 85 } 85 } 86 /* 86 /* 87 ********************************************** 87 ************************************************************ 88 */ 88 */ 89 static int MCGIDI_energyAngular_linear_parseFr 89 static int MCGIDI_energyAngular_linear_parseFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution ) { 90 90 91 int iV, iW; 91 int iV, iW; 92 double y, norm, energyInFactor, energyOutF 92 double y, norm, energyInFactor, energyOutFactor; 93 char const *energyUnit, *multiplicityProba 93 char const *energyUnit, *multiplicityProbabilityUnits[2] = { "", "1/MeV" }; 94 xDataTOM_element *linear; 94 xDataTOM_element *linear; 95 ptwXY_interpolation interpolationXY, inter 95 ptwXY_interpolation interpolationXY, interpolationWY, interpolationVY; 96 xDataTOM_XYs *XYs; 96 xDataTOM_XYs *XYs; 97 xDataTOM_W_XYs *W_XYs; 97 xDataTOM_W_XYs *W_XYs; 98 xDataTOM_V_W_XYs *V_W_XYs; 98 xDataTOM_V_W_XYs *V_W_XYs; 99 MCGIDI_pdfsOfXGivenW *pdfOfEpGivenE, *pdfO 99 MCGIDI_pdfsOfXGivenW *pdfOfEpGivenE, *pdfOfMuGivenEAndEp = NULL, *pdfOfMuGivenEAndEp2 = NULL; 100 MCGIDI_energyAngular *energyAngular = NULL 100 MCGIDI_energyAngular *energyAngular = NULL; 101 ptwXYPoints *pdfXY1 = NULL, *pdfXY2 = NULL 101 ptwXYPoints *pdfXY1 = NULL, *pdfXY2 = NULL; 102 nfu_status status; 102 nfu_status status; 103 103 104 if( ( linear = xDataTOME_getOneElementByNa 104 if( ( linear = xDataTOME_getOneElementByName( smr, element, "linear", 1 ) ) == NULL ) goto err; 105 105 106 if( MCGIDI_fromTOM_interpolation( smr, lin 106 if( MCGIDI_fromTOM_interpolation( smr, linear, 0, &interpolationVY ) ) goto err; 107 if( MCGIDI_fromTOM_interpolation( smr, lin 107 if( MCGIDI_fromTOM_interpolation( smr, linear, 1, &interpolationWY ) ) goto err; 108 if( MCGIDI_fromTOM_interpolation( smr, lin 108 if( MCGIDI_fromTOM_interpolation( smr, linear, 2, &interpolationXY ) ) goto err; 109 109 110 if( ( energyAngular = MCGIDI_energyAngular 110 if( ( energyAngular = MCGIDI_energyAngular_new( smr ) ) == NULL ) goto err; 111 if( ( energyAngular->frame = MCGIDI_misc_g 111 if( ( energyAngular->frame = MCGIDI_misc_getProductFrame( smr, linear ) ) == xDataTOM_frame_invalid ) goto err; 112 112 113 pdfOfEpGivenE = &(energyAngular->pdfOfEpGi 113 pdfOfEpGivenE = &(energyAngular->pdfOfEpGivenE); 114 pdfOfEpGivenE->interpolationWY = interpola 114 pdfOfEpGivenE->interpolationWY = interpolationVY; 115 pdfOfEpGivenE->interpolationXY = interpola 115 pdfOfEpGivenE->interpolationXY = interpolationWY; 116 116 117 if( ( V_W_XYs = (xDataTOM_V_W_XYs *) xData 117 if( ( V_W_XYs = (xDataTOM_V_W_XYs *) xDataTOME_getXDataIfID( smr, linear, "V_W_XYs" ) ) == NULL ) goto err; 118 if( ( pdfOfEpGivenE->Ws = (double *) smr_m 118 if( ( pdfOfEpGivenE->Ws = (double *) smr_malloc2( smr, V_W_XYs->length * sizeof( double ), 1, "pdfOfEpGivenE->Ws" ) ) == NULL ) goto err; 119 if( ( pdfOfEpGivenE->dist = (MCGIDI_pdfOfX 119 if( ( pdfOfEpGivenE->dist = (MCGIDI_pdfOfX *) smr_malloc2( smr, V_W_XYs->length * sizeof( MCGIDI_pdfOfX ), 0, "pdfOfEpGivenE->dist" ) ) == NULL ) goto err; 120 if( ( pdfOfMuGivenEAndEp = (MCGIDI_pdfsOfX 120 if( ( pdfOfMuGivenEAndEp = (MCGIDI_pdfsOfXGivenW *) smr_malloc2( smr, V_W_XYs->length * sizeof( MCGIDI_pdfsOfXGivenW ), 1, "pdfOfMuGivenEAndEp" ) ) == NULL ) goto err; 121 121 122 energyUnit = xDataTOM_subAxes_getUnit( smr 122 energyUnit = xDataTOM_subAxes_getUnit( smr, &(V_W_XYs->subAxes), 0 ); 123 if( !smr_isOk( smr ) ) goto err; 123 if( !smr_isOk( smr ) ) goto err; 124 energyInFactor = MCGIDI_misc_getUnitConver 124 energyInFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyUnit, "MeV" ); 125 if( !smr_isOk( smr ) ) goto err; 125 if( !smr_isOk( smr ) ) goto err; 126 126 127 energyUnit = xDataTOM_subAxes_getUnit( smr 127 energyUnit = xDataTOM_subAxes_getUnit( smr, &(V_W_XYs->subAxes), 1 ); 128 if( !smr_isOk( smr ) ) goto err; 128 if( !smr_isOk( smr ) ) goto err; 129 energyOutFactor = MCGIDI_misc_getUnitConve 129 energyOutFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyUnit, "MeV" ); 130 if( !smr_isOk( smr ) ) goto err; 130 if( !smr_isOk( smr ) ) goto err; 131 131 132 for( iV = 0; iV < V_W_XYs->length; iV++ ) 132 for( iV = 0; iV < V_W_XYs->length; iV++ ) { 133 W_XYs = &(V_W_XYs->W_XYs[iV]); 133 W_XYs = &(V_W_XYs->W_XYs[iV]); 134 pdfOfMuGivenEAndEp2 = &(pdfOfMuGivenEA 134 pdfOfMuGivenEAndEp2 = &(pdfOfMuGivenEAndEp[iV]); 135 pdfOfMuGivenEAndEp2->interpolationWY = 135 pdfOfMuGivenEAndEp2->interpolationWY = interpolationWY; 136 pdfOfMuGivenEAndEp2->interpolationXY = 136 pdfOfMuGivenEAndEp2->interpolationXY = interpolationXY; 137 if( ( pdfXY2 = ptwXY_new( interpolatio 137 if( ( pdfXY2 = ptwXY_new( interpolationWY, NULL, 2., 1e-6, W_XYs->length, 10, &status, 0 ) ) == NULL ) goto errA; 138 if( ( pdfOfMuGivenEAndEp2->Ws = (doubl 138 if( ( pdfOfMuGivenEAndEp2->Ws = (double *) smr_malloc2( smr, W_XYs->length * sizeof( double ), 1, "pdfOfMuGivenEAndEp2->Ws" ) ) == NULL ) goto err; 139 if( ( pdfOfMuGivenEAndEp2->dist = (MCG 139 if( ( pdfOfMuGivenEAndEp2->dist = (MCGIDI_pdfOfX *) smr_malloc2( smr, W_XYs->length * sizeof( MCGIDI_pdfOfX ), 0, "pdfOfMuGivenEAndEp2->dist" ) ) == NULL ) goto err; 140 for( iW = 0; iW < W_XYs->length; iW++ 140 for( iW = 0; iW < W_XYs->length; iW++ ) { 141 XYs = &(W_XYs->XYs[iW]); 141 XYs = &(W_XYs->XYs[iW]); 142 if( ( pdfXY1 = MCGIDI_misc_dataFr 142 if( ( pdfXY1 = MCGIDI_misc_dataFromXYs2ptwXYPointsInUnitsOf( smr, XYs, interpolationXY, multiplicityProbabilityUnits ) ) == NULL ) goto err; 143 y = ptwXY_integrateDomain( pdfXY1, 143 y = ptwXY_integrateDomain( pdfXY1, &status ); 144 if( ( status = ptwXY_setValueAtX( 144 if( ( status = ptwXY_setValueAtX( pdfXY2, energyOutFactor * XYs->value, y ) ) != nfu_Okay ) goto errA; >> 145 if( status != nfu_Okay ) goto errA; 145 146 146 if( y == 0 ) { 147 if( y == 0 ) { 147 if( ( status = ptwXY_add_doubl 148 if( ( status = ptwXY_add_double( pdfXY1, 0.5 ) ) != nfu_Okay ) goto errA; 148 } 149 } 149 pdfOfMuGivenEAndEp2->Ws[iW] = ener 150 pdfOfMuGivenEAndEp2->Ws[iW] = energyOutFactor * XYs->value; 150 if( MCGIDI_fromTOM_pdfOfX( smr, pd 151 if( MCGIDI_fromTOM_pdfOfX( smr, pdfXY1, &(pdfOfMuGivenEAndEp2->dist[iW]), &norm ) ) goto err; 151 pdfOfMuGivenEAndEp2->numberOfWs++; 152 pdfOfMuGivenEAndEp2->numberOfWs++; 152 153 153 pdfXY1 = ptwXY_free( pdfXY1 ); 154 pdfXY1 = ptwXY_free( pdfXY1 ); 154 } 155 } 155 pdfOfEpGivenE->Ws[iV] = energyInFactor 156 pdfOfEpGivenE->Ws[iV] = energyInFactor * W_XYs->value; 156 if( MCGIDI_fromTOM_pdfOfX( smr, pdfXY2 157 if( MCGIDI_fromTOM_pdfOfX( smr, pdfXY2, &(pdfOfEpGivenE->dist[iV]), &norm ) ) goto err; 157 pdfOfEpGivenE->numberOfWs++; 158 pdfOfEpGivenE->numberOfWs++; 158 159 159 pdfXY2 = ptwXY_free( pdfXY2 ); 160 pdfXY2 = ptwXY_free( pdfXY2 ); 160 } 161 } 161 energyAngular->pdfOfMuGivenEAndEp = pdfOfM 162 energyAngular->pdfOfMuGivenEAndEp = pdfOfMuGivenEAndEp; 162 distribution->energyAngular = energyAngula 163 distribution->energyAngular = energyAngular; 163 distribution->type = MCGIDI_distributionTy 164 distribution->type = MCGIDI_distributionType_energyAngular_e; 164 165 165 return( 0 ); 166 return( 0 ); 166 167 167 errA: 168 errA: 168 smr_setReportError2( smr, smr_unknownID, 1 169 smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_integrateDomain err = %d: %s\n", status, nfu_statusMessage( status ) ); 169 err: 170 err: 170 if( pdfXY1 != NULL ) ptwXY_free( pdfXY1 ); 171 if( pdfXY1 != NULL ) ptwXY_free( pdfXY1 ); 171 if( pdfXY2 != NULL ) ptwXY_free( pdfXY2 ); 172 if( pdfXY2 != NULL ) ptwXY_free( pdfXY2 ); 172 if( energyAngular != NULL ) MCGIDI_energyA 173 if( energyAngular != NULL ) MCGIDI_energyAngular_free( smr, energyAngular ); 173 /* ????????? Need to free pdfOfMuGivenEAndEp, 174 /* ????????? Need to free pdfOfMuGivenEAndEp, now may be handled by MCGIDI_energyAngular_free. Need to check. */ 174 return( 1 ); 175 return( 1 ); 175 } 176 } 176 /* 177 /* 177 ********************************************** 178 ************************************************************ 178 */ 179 */ 179 int MCGIDI_energyAngular_sampleDistribution( s 180 int MCGIDI_energyAngular_sampleDistribution( statusMessageReporting *smr, MCGIDI_distribution *distribution, MCGIDI_quantitiesLookupModes &modes, 180 MCGIDI_decaySamplingInfo *decaySamplin 181 MCGIDI_decaySamplingInfo *decaySamplingInfo ) { 181 182 182 double Ep; 183 double Ep; 183 MCGIDI_energyAngular *energyAngular = dist 184 MCGIDI_energyAngular *energyAngular = distribution->energyAngular; 184 185 185 MCGIDI_sampling_doubleDistribution( smr, & 186 MCGIDI_sampling_doubleDistribution( smr, &(energyAngular->pdfOfEpGivenE), energyAngular->pdfOfMuGivenEAndEp, modes, decaySamplingInfo ); 186 Ep = decaySamplingInfo->mu; 187 Ep = decaySamplingInfo->mu; 187 decaySamplingInfo->mu = decaySamplingInfo- 188 decaySamplingInfo->mu = decaySamplingInfo->Ep; 188 decaySamplingInfo->Ep = Ep; 189 decaySamplingInfo->Ep = Ep; 189 decaySamplingInfo->frame = energyAngular-> 190 decaySamplingInfo->frame = energyAngular->frame; 190 191 191 return( 0 ); 192 return( 0 ); 192 } 193 } 193 194 194 #if defined __cplusplus 195 #if defined __cplusplus 195 } 196 } 196 #endif 197 #endif 197 198 198 199