Geant4 Cross Reference |
1 /* 1 /* 2 # <<BEGIN-copyright>> 2 # <<BEGIN-copyright>> 3 # <<END-copyright>> 3 # <<END-copyright>> 4 */ 4 */ 5 #include <stdio.h> 5 #include <stdio.h> 6 #include <stdlib.h> 6 #include <stdlib.h> 7 #include <string.h> 7 #include <string.h> 8 #include <cmath> 8 #include <cmath> 9 /* 9 /* 10 #include <unistd.h> 10 #include <unistd.h> 11 #include <ctype.h> 11 #include <ctype.h> 12 */ 12 */ 13 13 14 #include "MCGIDI_fromTOM.h" 14 #include "MCGIDI_fromTOM.h" 15 #include "MCGIDI_misc.h" 15 #include "MCGIDI_misc.h" 16 16 17 #if defined __cplusplus 17 #if defined __cplusplus 18 namespace GIDI { 18 namespace GIDI { 19 using namespace GIDI; 19 using namespace GIDI; 20 #endif 20 #endif 21 21 22 static int MCGIDI_fromTOM_pdfOfXGivenW( status 22 static int MCGIDI_fromTOM_pdfOfXGivenW( statusMessageReporting *smr, ptwXYPoints *pdfXY, MCGIDI_pdfsOfXGivenW *dists, int i, double *norm ); 23 /* 23 /* 24 ********************************************** 24 ************************************************************ 25 */ 25 */ 26 int MCGIDI_fromTOM_pdfsOfXGivenW( statusMessag 26 int MCGIDI_fromTOM_pdfsOfXGivenW( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_pdfsOfXGivenW *dists, ptwXYPoints *norms, 27 char const *toUnits[3] ) { 27 char const *toUnits[3] ) { 28 28 29 int i; 29 int i; 30 double norm, wUnitFactor; 30 double norm, wUnitFactor; 31 char const *wFromUnit, *toUnitsXY[2] = { t 31 char const *wFromUnit, *toUnitsXY[2] = { toUnits[1], toUnits[2] }; 32 xDataTOM_XYs *XYs; 32 xDataTOM_XYs *XYs; 33 xDataTOM_W_XYs *W_XYs; 33 xDataTOM_W_XYs *W_XYs; 34 ptwXYPoints *pdfXY = NULL; 34 ptwXYPoints *pdfXY = NULL; 35 ptwXY_interpolation interpolationXY, inter 35 ptwXY_interpolation interpolationXY, interpolationWY; 36 36 37 wFromUnit = xDataTOM_axes_getUnit( smr, &( 37 wFromUnit = xDataTOM_axes_getUnit( smr, &(element->xDataInfo.axes), 0 ); 38 if( !smr_isOk( smr ) ) goto err; 38 if( !smr_isOk( smr ) ) goto err; 39 wUnitFactor = MCGIDI_misc_getUnitConversio 39 wUnitFactor = MCGIDI_misc_getUnitConversionFactor( smr, wFromUnit, toUnits[0] ); 40 if( !smr_isOk( smr ) ) goto err; 40 if( !smr_isOk( smr ) ) goto err; 41 41 42 if( MCGIDI_fromTOM_interpolation( smr, ele 42 if( MCGIDI_fromTOM_interpolation( smr, element, 0, &interpolationWY ) ) goto err; 43 if( MCGIDI_fromTOM_interpolation( smr, ele 43 if( MCGIDI_fromTOM_interpolation( smr, element, 1, &interpolationXY ) ) goto err; 44 dists->interpolationWY = interpolationWY; 44 dists->interpolationWY = interpolationWY; 45 dists->interpolationXY = interpolationXY; 45 dists->interpolationXY = interpolationXY; 46 if( norms != NULL ) { 46 if( norms != NULL ) { 47 if( interpolationWY == ptwXY_interpola 47 if( interpolationWY == ptwXY_interpolationOther ) { 48 smr_setReportError2p( smr, smr_unk 48 smr_setReportError2p( smr, smr_unknownID, 1, "interpolationWY ptwXY_interpolationOther not supported" ); 49 goto err; 49 goto err; 50 } 50 } 51 } 51 } 52 52 53 W_XYs = (xDataTOM_W_XYs *) xDataTOME_getXD 53 W_XYs = (xDataTOM_W_XYs *) xDataTOME_getXDataIfID( smr, element, "W_XYs" ); 54 if( ( dists->Ws = (double *) smr_malloc2( 54 if( ( dists->Ws = (double *) smr_malloc2( smr, W_XYs->length * sizeof( double ), 1, "dists->Ws" ) ) == NULL ) goto err; 55 if( ( dists->dist = (MCGIDI_pdfOfX *) smr_ 55 if( ( dists->dist = (MCGIDI_pdfOfX *) smr_malloc2( smr, W_XYs->length * sizeof( MCGIDI_pdfOfX ), 0, "dists->dist" ) ) == NULL ) goto err; 56 56 57 for( i = 0; i < W_XYs->length; i++ ) { 57 for( i = 0; i < W_XYs->length; i++ ) { 58 XYs = &(W_XYs->XYs[i]); 58 XYs = &(W_XYs->XYs[i]); 59 dists->Ws[i] = wUnitFactor * XYs->valu 59 dists->Ws[i] = wUnitFactor * XYs->value; 60 if( ( pdfXY = MCGIDI_misc_dataFromXYs 60 if( ( pdfXY = MCGIDI_misc_dataFromXYs2ptwXYPointsInUnitsOf( smr, XYs, interpolationXY, toUnitsXY ) ) == NULL ) goto err; 61 if( MCGIDI_fromTOM_pdfOfXGivenW( smr, 61 if( MCGIDI_fromTOM_pdfOfXGivenW( smr, pdfXY, dists, i, &norm ) ) goto err; 62 if( norms != NULL ) { 62 if( norms != NULL ) { 63 ptwXY_setValueAtX( norms, XYs->val 63 ptwXY_setValueAtX( norms, XYs->value, norm ); } 64 else if( std::fabs( 1. - norm ) > 0.99 64 else if( std::fabs( 1. - norm ) > 0.99 ) { 65 smr_setReportError2( smr, smr_unkn 65 smr_setReportError2( smr, smr_unknownID, 1, "bad norm = %e for data", norm ); 66 goto err; 66 goto err; 67 } 67 } 68 ptwXY_free( pdfXY ); 68 ptwXY_free( pdfXY ); 69 pdfXY = NULL; 69 pdfXY = NULL; 70 } 70 } 71 71 72 return( 0 ); 72 return( 0 ); 73 73 74 err: 74 err: 75 if( pdfXY != NULL ) ptwXY_free( pdfXY ); 75 if( pdfXY != NULL ) ptwXY_free( pdfXY ); 76 return( 1 ); 76 return( 1 ); 77 } 77 } 78 /* 78 /* 79 ********************************************** 79 ************************************************************ 80 */ 80 */ 81 static int MCGIDI_fromTOM_pdfOfXGivenW( status 81 static int MCGIDI_fromTOM_pdfOfXGivenW( statusMessageReporting *smr, ptwXYPoints *pdfXY, MCGIDI_pdfsOfXGivenW *dists, int i, double *norm ) { 82 82 83 if( MCGIDI_fromTOM_pdfOfX( smr, pdfXY, &(d 83 if( MCGIDI_fromTOM_pdfOfX( smr, pdfXY, &(dists->dist[i]), norm ) ) return( 1 ); 84 dists->numberOfWs++; 84 dists->numberOfWs++; 85 return( 0 ); 85 return( 0 ); 86 } 86 } 87 /* 87 /* 88 ********************************************** 88 ************************************************************ 89 */ 89 */ 90 int MCGIDI_fromTOM_pdfOfX( statusMessageReport 90 int MCGIDI_fromTOM_pdfOfX( statusMessageReporting *smr, ptwXYPoints *pdfXY, MCGIDI_pdfOfX *dist, double *norm ) { 91 91 92 int j1, n1 = (int) ptwXY_length( pdfXY ); 92 int j1, n1 = (int) ptwXY_length( pdfXY ); 93 nfu_status status; 93 nfu_status status; 94 ptwXPoints *cdfX = NULL; 94 ptwXPoints *cdfX = NULL; 95 ptwXYPoint *point; 95 ptwXYPoint *point; 96 96 97 dist->numberOfXs = 0; 97 dist->numberOfXs = 0; 98 dist->Xs = NULL; 98 dist->Xs = NULL; 99 if( ptwXY_simpleCoalescePoints( pdfXY ) != 99 if( ptwXY_simpleCoalescePoints( pdfXY ) != nfu_Okay ) goto err; 100 100 101 if( ( dist->Xs = (double *) smr_malloc2( s 101 if( ( dist->Xs = (double *) smr_malloc2( smr, 3 * n1 * sizeof( double ), 0, "dist->Xs" ) ) == NULL ) goto err; 102 dist->pdf = &(dist->Xs[n1]); 102 dist->pdf = &(dist->Xs[n1]); 103 dist->cdf = &(dist->pdf[n1]); 103 dist->cdf = &(dist->pdf[n1]); 104 104 105 for( j1 = 0; j1 < n1; j1++ ) { 105 for( j1 = 0; j1 < n1; j1++ ) { 106 point = ptwXY_getPointAtIndex_Unsa 106 point = ptwXY_getPointAtIndex_Unsafely( pdfXY, j1 ); 107 dist->Xs[j1] = point->x; 107 dist->Xs[j1] = point->x; 108 dist->pdf[j1] = point->y; 108 dist->pdf[j1] = point->y; 109 } 109 } 110 110 111 if( ( cdfX = ptwXY_runningIntegral( pdfXY, 111 if( ( cdfX = ptwXY_runningIntegral( pdfXY, &status ) ) == NULL ) { 112 smr_setReportError2( smr, smr_unkn 112 smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_runningIntegral err = %d: %s\n", status, nfu_statusMessage( status ) ); 113 goto err; 113 goto err; 114 } 114 } 115 *norm = ptwX_getPointAtIndex_Unsafely( cdf 115 *norm = ptwX_getPointAtIndex_Unsafely( cdfX, n1 - 1 ); 116 if( *norm == 0. ) { /* Should 116 if( *norm == 0. ) { /* Should only happend for gammas. */ 117 double inv_norm, sum = 0; 117 double inv_norm, sum = 0; 118 118 119 inv_norm = 1.0 / ( dist->Xs[n1-1] - di 119 inv_norm = 1.0 / ( dist->Xs[n1-1] - dist->Xs[0] ); 120 for( j1 = 0; j1 < n1; ++j1 ) { 120 for( j1 = 0; j1 < n1; ++j1 ) { 121 if( j1 > 0 ) sum += dist->Xs[j1] - 121 if( j1 > 0 ) sum += dist->Xs[j1] - dist->Xs[j1-1]; 122 dist->pdf[j1] = 1; 122 dist->pdf[j1] = 1; 123 dist->cdf[j1] = sum * inv_norm; 123 dist->cdf[j1] = sum * inv_norm; 124 } 124 } 125 dist->cdf[n1-1] = 1.; } 125 dist->cdf[n1-1] = 1.; } 126 else { 126 else { 127 for( j1 = 0; j1 < n1; j1++ ) dist->cdf 127 for( j1 = 0; j1 < n1; j1++ ) dist->cdf[j1] = ptwX_getPointAtIndex_Unsafely( cdfX, j1 ) / *norm; 128 for( j1 = 0; j1 < n1; j1++ ) dist->pdf 128 for( j1 = 0; j1 < n1; j1++ ) dist->pdf[j1] /= *norm; 129 } 129 } 130 ptwX_free( cdfX ); 130 ptwX_free( cdfX ); 131 131 132 dist->numberOfXs = n1; 132 dist->numberOfXs = n1; 133 return( 0 ); 133 return( 0 ); 134 134 135 err: 135 err: 136 if( dist->Xs != NULL ) smr_freeMemory( (vo 136 if( dist->Xs != NULL ) smr_freeMemory( (void **) &(dist->Xs) ); 137 if( cdfX != NULL ) ptwX_free( cdfX ); 137 if( cdfX != NULL ) ptwX_free( cdfX ); 138 return( 1 ); 138 return( 1 ); 139 } 139 } 140 /* 140 /* 141 ********************************************** 141 ************************************************************ 142 */ 142 */ 143 int MCGIDI_fromTOM_interpolation( statusMessag 143 int MCGIDI_fromTOM_interpolation( statusMessageReporting *smr, xDataTOM_element *element, int index, ptwXY_interpolation *interpolation ) { 144 144 145 enum xDataTOM_interpolationFlag independen 145 enum xDataTOM_interpolationFlag independent, dependent; 146 enum xDataTOM_interpolationQualifier quali 146 enum xDataTOM_interpolationQualifier qualifier; 147 147 148 if( xDataTOME_getInterpolation( smr, eleme 148 if( xDataTOME_getInterpolation( smr, element, index, &independent, &dependent, &qualifier ) ) return( 1 ); 149 149 150 *interpolation = ptwXY_interpolationOther; 150 *interpolation = ptwXY_interpolationOther; 151 151 152 if( dependent == xDataTOM_interpolationFla 152 if( dependent == xDataTOM_interpolationFlag_flat ) { 153 *interpolation = ptwXY_interpolationFl 153 *interpolation = ptwXY_interpolationFlat; } 154 else if( independent == xDataTOM_interpola 154 else if( independent == xDataTOM_interpolationFlag_linear ) { 155 if( dependent == xDataTOM_interpolatio 155 if( dependent == xDataTOM_interpolationFlag_linear ) { 156 *interpolation = ptwXY_interpolati 156 *interpolation = ptwXY_interpolationLinLin; } 157 else if( dependent == xDataTOM_interpo 157 else if( dependent == xDataTOM_interpolationFlag_log ) { 158 *interpolation = ptwXY_interpolati 158 *interpolation = ptwXY_interpolationLinLog; 159 } } 159 } } 160 else if( independent == xDataTOM_interpola 160 else if( independent == xDataTOM_interpolationFlag_log ) { 161 if( dependent == xDataTOM_interpolatio 161 if( dependent == xDataTOM_interpolationFlag_linear ) { 162 *interpolation = ptwXY_interpolati 162 *interpolation = ptwXY_interpolationLogLin; } 163 else if( dependent == xDataTOM_interpo 163 else if( dependent == xDataTOM_interpolationFlag_log ) { 164 *interpolation = ptwXY_interpolati 164 *interpolation = ptwXY_interpolationLogLog; 165 } 165 } 166 } 166 } 167 167 168 return( 0 ); 168 return( 0 ); 169 } 169 } 170 170 171 #if defined __cplusplus 171 #if defined __cplusplus 172 } 172 } 173 #endif 173 #endif 174 174