Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lend/src/MCGIDI_fromTOM.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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