Geant4 Cross Reference

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