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.h" 7 #include "MCGIDI.h" 8 #include "MCGIDI_fromTOM.h" 8 #include "MCGIDI_fromTOM.h" 9 #include "MCGIDI_misc.h" 9 #include "MCGIDI_misc.h" 10 #include "MCGIDI_private.h" 10 #include "MCGIDI_private.h" 11 11 12 #if defined __cplusplus 12 #if defined __cplusplus 13 namespace GIDI { 13 namespace GIDI { 14 using namespace GIDI; 14 using namespace GIDI; 15 #endif 15 #endif 16 16 17 static int MCGIDI_angularEnergy_parsePointwise 17 static int MCGIDI_angularEnergy_parsePointwiseFromTOM( statusMessageReporting *smr, xDataTOM_element *pointwise, MCGIDI_distribution *distribution ); 18 /* 18 /* 19 ********************************************** 19 ************************************************************ 20 */ 20 */ 21 MCGIDI_angularEnergy *MCGIDI_angularEnergy_new 21 MCGIDI_angularEnergy *MCGIDI_angularEnergy_new( statusMessageReporting *smr ) { 22 22 23 MCGIDI_angularEnergy *angularEnergy; 23 MCGIDI_angularEnergy *angularEnergy; 24 24 25 if( ( angularEnergy = (MCGIDI_angularEnerg 25 if( ( angularEnergy = (MCGIDI_angularEnergy *) smr_malloc2( smr, sizeof( MCGIDI_angularEnergy ), 0, "angularEnergy" ) ) == NULL ) return( NULL ); 26 if( MCGIDI_angularEnergy_initialize( smr, 26 if( MCGIDI_angularEnergy_initialize( smr, angularEnergy ) ) angularEnergy = MCGIDI_angularEnergy_free( smr, angularEnergy ); 27 return( angularEnergy ); 27 return( angularEnergy ); 28 } 28 } 29 /* 29 /* 30 ********************************************** 30 ************************************************************ 31 */ 31 */ 32 int MCGIDI_angularEnergy_initialize( statusMes 32 int MCGIDI_angularEnergy_initialize( statusMessageReporting * /*smr*/, MCGIDI_angularEnergy *angularEnergy ) { 33 33 34 memset( angularEnergy, 0, sizeof( MCGIDI_a 34 memset( angularEnergy, 0, sizeof( MCGIDI_angularEnergy ) ); 35 return( 0 ); 35 return( 0 ); 36 } 36 } 37 /* 37 /* 38 ********************************************** 38 ************************************************************ 39 */ 39 */ 40 MCGIDI_angularEnergy *MCGIDI_angularEnergy_fre 40 MCGIDI_angularEnergy *MCGIDI_angularEnergy_free( statusMessageReporting *smr, MCGIDI_angularEnergy *angularEnergy ) { 41 41 42 MCGIDI_angularEnergy_release( smr, angular 42 MCGIDI_angularEnergy_release( smr, angularEnergy ); 43 smr_freeMemory( (void **) &angularEnergy ) 43 smr_freeMemory( (void **) &angularEnergy ); 44 return( NULL ); 44 return( NULL ); 45 } 45 } 46 /* 46 /* 47 ********************************************** 47 ************************************************************ 48 */ 48 */ 49 int MCGIDI_angularEnergy_release( statusMessag 49 int MCGIDI_angularEnergy_release( statusMessageReporting *smr, MCGIDI_angularEnergy *angularEnergy ) { 50 50 51 int i; 51 int i; 52 52 53 for( i = 0; i < angularEnergy->pdfOfMuGive 53 for( i = 0; i < angularEnergy->pdfOfMuGivenE.numberOfWs; i++ ) { 54 MCGIDI_sampling_pdfsOfXGivenW_release( 54 MCGIDI_sampling_pdfsOfXGivenW_release( smr, &(angularEnergy->pdfOfEpGivenEAndMu[i]) ); 55 } 55 } 56 smr_freeMemory( (void **) &(angularEnergy- 56 smr_freeMemory( (void **) &(angularEnergy->pdfOfEpGivenEAndMu) ); 57 MCGIDI_sampling_pdfsOfXGivenW_release( smr 57 MCGIDI_sampling_pdfsOfXGivenW_release( smr, &(angularEnergy->pdfOfMuGivenE) ); 58 MCGIDI_angularEnergy_initialize( smr, angu 58 MCGIDI_angularEnergy_initialize( smr, angularEnergy ); 59 59 60 return( 0 ); 60 return( 0 ); 61 } 61 } 62 /* 62 /* 63 ********************************************** 63 ************************************************************ 64 */ 64 */ 65 int MCGIDI_angularEnergy_parseFromTOM( statusM 65 int MCGIDI_angularEnergy_parseFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution ) { 66 66 67 xDataTOM_element *angularEnergyElement, *p 67 xDataTOM_element *angularEnergyElement, *pointwise = NULL; 68 char const *nativeData; 68 char const *nativeData; 69 69 70 if( ( angularEnergyElement = xDataTOME_get 70 if( ( angularEnergyElement = xDataTOME_getOneElementByName( smr, element, "angularEnergy", 1 ) ) == NULL ) goto err; 71 71 72 if( ( nativeData = xDataTOM_getAttributesV 72 if( ( nativeData = xDataTOM_getAttributesValueInElement( angularEnergyElement, "nativeData" ) ) == NULL ) goto err; 73 if( strcmp( nativeData, "pointwise" ) == 0 73 if( strcmp( nativeData, "pointwise" ) == 0 ) { 74 if( ( pointwise = xDataTOME_getOneElem 74 if( ( pointwise = xDataTOME_getOneElementByName( smr, angularEnergyElement, "pointwise", 1 ) ) == NULL ) goto err; } 75 else if( strcmp( nativeData, "linear" ) == 75 else if( strcmp( nativeData, "linear" ) == 0 ) { 76 if( ( pointwise = xDataTOME_getOneElem 76 if( ( pointwise = xDataTOME_getOneElementByName( smr, angularEnergyElement, "linear", 1 ) ) == NULL ) goto err; } 77 else { 77 else { 78 smr_setReportError2( smr, smr_unknownI 78 smr_setReportError2( smr, smr_unknownID, 1, "angularEnergy nativeData = '%s' not supported", nativeData ); 79 goto err; 79 goto err; 80 } 80 } 81 if( pointwise != NULL ) return( MCGIDI_ang 81 if( pointwise != NULL ) return( MCGIDI_angularEnergy_parsePointwiseFromTOM( smr, pointwise, distribution ) ); 82 82 83 return( 0 ); 83 return( 0 ); 84 84 85 err: 85 err: 86 return( 1 ); 86 return( 1 ); 87 } 87 } 88 /* 88 /* 89 ********************************************** 89 ************************************************************ 90 */ 90 */ 91 static int MCGIDI_angularEnergy_parsePointwise 91 static int MCGIDI_angularEnergy_parsePointwiseFromTOM( statusMessageReporting *smr, xDataTOM_element *pointwise, MCGIDI_distribution *distribution ) { 92 92 93 int iV, iW; 93 int iV, iW; 94 double y, norm, energyInFactor; 94 double y, norm, energyInFactor; 95 char const *energyUnit, *energyOutProbabil 95 char const *energyUnit, *energyOutProbabilityUnits[2] = { "MeV", "1/MeV" }; 96 MCGIDI_angularEnergy *angularEnergy = NULL 96 MCGIDI_angularEnergy *angularEnergy = NULL; 97 ptwXY_interpolation interpolationXY, inter 97 ptwXY_interpolation interpolationXY, interpolationWY, interpolationVY; 98 xDataTOM_XYs *XYs; 98 xDataTOM_XYs *XYs; 99 xDataTOM_W_XYs *W_XYs; 99 xDataTOM_W_XYs *W_XYs; 100 xDataTOM_V_W_XYs *V_W_XYs; 100 xDataTOM_V_W_XYs *V_W_XYs; 101 MCGIDI_pdfsOfXGivenW *pdfOfMuGivenE, *pdfO 101 MCGIDI_pdfsOfXGivenW *pdfOfMuGivenE, *pdfOfEpGivenEAndMu = NULL, *pdfOfEpGivenEAndMu2 = NULL; 102 ptwXYPoints *pdfXY1 = NULL, *pdfXY2 = NULL 102 ptwXYPoints *pdfXY1 = NULL, *pdfXY2 = NULL; 103 nfu_status status; 103 nfu_status status; 104 104 105 if( MCGIDI_fromTOM_interpolation( smr, poi 105 if( MCGIDI_fromTOM_interpolation( smr, pointwise, 0, &interpolationVY ) ) goto err; 106 if( MCGIDI_fromTOM_interpolation( smr, poi 106 if( MCGIDI_fromTOM_interpolation( smr, pointwise, 1, &interpolationWY ) ) goto err; 107 if( MCGIDI_fromTOM_interpolation( smr, poi 107 if( MCGIDI_fromTOM_interpolation( smr, pointwise, 2, &interpolationXY ) ) goto err; 108 if( ( angularEnergy = MCGIDI_angularEnergy 108 if( ( angularEnergy = MCGIDI_angularEnergy_new( smr ) ) == NULL ) goto err; 109 109 110 if( ( angularEnergy->frame = MCGIDI_misc_g 110 if( ( angularEnergy->frame = MCGIDI_misc_getProductFrame( smr, pointwise ) ) == xDataTOM_frame_invalid ) goto err; 111 111 112 pdfOfMuGivenE = &(angularEnergy->pdfOfMuGi 112 pdfOfMuGivenE = &(angularEnergy->pdfOfMuGivenE); 113 pdfOfMuGivenE->interpolationWY = interpola 113 pdfOfMuGivenE->interpolationWY = interpolationVY; 114 pdfOfMuGivenE->interpolationXY = interpola 114 pdfOfMuGivenE->interpolationXY = interpolationWY; 115 115 116 if( ( V_W_XYs = (xDataTOM_V_W_XYs *) xData 116 if( ( V_W_XYs = (xDataTOM_V_W_XYs *) xDataTOME_getXDataIfID( smr, pointwise, "V_W_XYs" ) ) == NULL ) goto err; 117 if( ( pdfOfMuGivenE->Ws = (double *) smr_m 117 if( ( pdfOfMuGivenE->Ws = (double *) smr_malloc2( smr, V_W_XYs->length * sizeof( double ), 1, "pdfOfMuGivenE->Ws" ) ) == NULL ) goto err; 118 if( ( pdfOfMuGivenE->dist = (MCGIDI_pdfOfX 118 if( ( pdfOfMuGivenE->dist = (MCGIDI_pdfOfX *) smr_malloc2( smr, V_W_XYs->length * sizeof( MCGIDI_pdfOfX ), 0, "pdfOfMuGivenE->dist" ) ) == NULL ) goto err; 119 if( ( pdfOfEpGivenEAndMu = (MCGIDI_pdfsOfX 119 if( ( pdfOfEpGivenEAndMu = (MCGIDI_pdfsOfXGivenW *) smr_malloc2( smr, V_W_XYs->length * sizeof( MCGIDI_pdfsOfXGivenW ), 1, "pdfOfEpGivenEAndMu" ) ) == NULL ) goto err; 120 120 121 energyUnit = xDataTOM_subAxes_getUnit( smr 121 energyUnit = xDataTOM_subAxes_getUnit( smr, &(V_W_XYs->subAxes), 0 ); 122 if( !smr_isOk( smr ) ) goto err; 122 if( !smr_isOk( smr ) ) goto err; 123 energyInFactor = MCGIDI_misc_getUnitConver 123 energyInFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyUnit, "MeV" ); 124 if( !smr_isOk( smr ) ) goto err; 124 if( !smr_isOk( smr ) ) goto err; 125 125 126 for( iV = 0; iV < V_W_XYs->length; iV++ ) 126 for( iV = 0; iV < V_W_XYs->length; iV++ ) { 127 W_XYs = &(V_W_XYs->W_XYs[iV]); 127 W_XYs = &(V_W_XYs->W_XYs[iV]); 128 pdfOfEpGivenEAndMu2 = &(pdfOfEpGivenEA 128 pdfOfEpGivenEAndMu2 = &(pdfOfEpGivenEAndMu[iV]); 129 pdfOfEpGivenEAndMu2->interpolationWY = 129 pdfOfEpGivenEAndMu2->interpolationWY = interpolationWY; 130 pdfOfEpGivenEAndMu2->interpolationXY = 130 pdfOfEpGivenEAndMu2->interpolationXY = interpolationXY; 131 if( ( pdfXY2 = ptwXY_new( interpolatio 131 if( ( pdfXY2 = ptwXY_new( interpolationWY, NULL, 2., 1e-6, W_XYs->length, 10, &status, 0 ) ) == NULL ) goto errA; 132 if( ( pdfOfEpGivenEAndMu2->Ws = (doubl 132 if( ( pdfOfEpGivenEAndMu2->Ws = (double *) smr_malloc2( smr, W_XYs->length * sizeof( double ), 1, "pdfOfEpGivenEAndMu2->Ws" ) ) == NULL ) goto err; 133 if( ( pdfOfEpGivenEAndMu2->dist = (MCG 133 if( ( pdfOfEpGivenEAndMu2->dist = (MCGIDI_pdfOfX *) smr_malloc2( smr, W_XYs->length * sizeof( MCGIDI_pdfOfX ), 0, "pdfOfEpGivenEAndMu2->dist" ) ) == NULL ) goto err; 134 for( iW = 0; iW < W_XYs->length; iW++ 134 for( iW = 0; iW < W_XYs->length; iW++ ) { 135 XYs = &(W_XYs->XYs[iW]); 135 XYs = &(W_XYs->XYs[iW]); 136 if( ( pdfXY1 = MCGIDI_misc_dataFr 136 if( ( pdfXY1 = MCGIDI_misc_dataFromXYs2ptwXYPointsInUnitsOf( smr, XYs, interpolationXY, energyOutProbabilityUnits ) ) == NULL ) goto err; 137 y = ptwXY_integrateDomain( pdfXY1, 137 y = ptwXY_integrateDomain( pdfXY1, &status ); 138 if( ( status = ptwXY_setValueAtX( 138 if( ( status = ptwXY_setValueAtX( pdfXY2, XYs->value, y ) ) != nfu_Okay ) goto errA; >> 139 if( status != nfu_Okay ) goto errA; 139 140 140 if( y == 0 ) { 141 if( y == 0 ) { 141 if( ( status = ptwXY_add_doubl 142 if( ( status = ptwXY_add_double( pdfXY1, 0.5 ) ) != nfu_Okay ) goto errA; 142 } 143 } 143 pdfOfEpGivenEAndMu2->Ws[iW] = XYs- 144 pdfOfEpGivenEAndMu2->Ws[iW] = XYs->value; 144 if( MCGIDI_fromTOM_pdfOfX( smr, pd 145 if( MCGIDI_fromTOM_pdfOfX( smr, pdfXY1, &(pdfOfEpGivenEAndMu2->dist[iW]), &norm ) ) goto err; 145 pdfOfEpGivenEAndMu2->numberOfWs++; 146 pdfOfEpGivenEAndMu2->numberOfWs++; 146 147 147 pdfXY1 = ptwXY_free( pdfXY1 ); 148 pdfXY1 = ptwXY_free( pdfXY1 ); 148 } 149 } 149 pdfOfMuGivenE->Ws[iV] = energyInFactor 150 pdfOfMuGivenE->Ws[iV] = energyInFactor * W_XYs->value; 150 if( MCGIDI_fromTOM_pdfOfX( smr, pdfXY2 151 if( MCGIDI_fromTOM_pdfOfX( smr, pdfXY2, &(pdfOfMuGivenE->dist[iV]), &norm ) ) goto err; 151 pdfOfMuGivenE->numberOfWs++; 152 pdfOfMuGivenE->numberOfWs++; 152 153 153 pdfXY2 = ptwXY_free( pdfXY2 ); 154 pdfXY2 = ptwXY_free( pdfXY2 ); 154 } 155 } 155 156 156 angularEnergy->pdfOfEpGivenEAndMu = pdfOfE 157 angularEnergy->pdfOfEpGivenEAndMu = pdfOfEpGivenEAndMu; 157 distribution->angularEnergy = angularEnerg 158 distribution->angularEnergy = angularEnergy; 158 distribution->type = MCGIDI_distributionTy 159 distribution->type = MCGIDI_distributionType_angularEnergy_e; 159 160 160 return( 0 ); 161 return( 0 ); 161 162 162 errA: 163 errA: 163 smr_setReportError2( smr, smr_unknownID, 1 164 smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_integrateDomain err = %d: %s\n", status, nfu_statusMessage( status ) ); 164 err: 165 err: 165 if( pdfXY1 != NULL ) ptwXY_free( pdfXY1 ); 166 if( pdfXY1 != NULL ) ptwXY_free( pdfXY1 ); 166 if( pdfXY2 != NULL ) ptwXY_free( pdfXY2 ); 167 if( pdfXY2 != NULL ) ptwXY_free( pdfXY2 ); 167 /* Need to free pdfOfEpGivenEAndMu. */ 168 /* Need to free pdfOfEpGivenEAndMu. */ 168 if( angularEnergy != NULL ) MCGIDI_angular 169 if( angularEnergy != NULL ) MCGIDI_angularEnergy_free( smr, angularEnergy ); 169 return( 1 ); 170 return( 1 ); 170 } 171 } 171 /* 172 /* 172 ********************************************** 173 ************************************************************ 173 */ 174 */ 174 int MCGIDI_angularEnergy_sampleDistribution( s 175 int MCGIDI_angularEnergy_sampleDistribution( statusMessageReporting *smr, MCGIDI_angularEnergy *angularEnergy, MCGIDI_quantitiesLookupModes &modes, 175 MCGIDI_decaySamplingInfo *decaySamplin 176 MCGIDI_decaySamplingInfo *decaySamplingInfo ) { 176 177 177 int status = MCGIDI_sampling_doubleDistrib 178 int status = MCGIDI_sampling_doubleDistribution( smr, &(angularEnergy->pdfOfMuGivenE), angularEnergy->pdfOfEpGivenEAndMu, modes, decaySamplingInfo ); 178 179 179 decaySamplingInfo->frame = angularEnergy-> 180 decaySamplingInfo->frame = angularEnergy->frame; 180 return( status ); 181 return( status ); 181 } 182 } 182 183 183 #if defined __cplusplus 184 #if defined __cplusplus 184 } 185 } 185 #endif 186 #endif 186 187 187 188