ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MCGIDI_angularEnergy.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MCGIDI_angularEnergy.cc
1 /*
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5 #include <string.h>
6 
7 #include "MCGIDI.h"
8 #include "MCGIDI_fromTOM.h"
9 #include "MCGIDI_misc.h"
10 #include "MCGIDI_private.h"
11 
12 #if defined __cplusplus
13 namespace GIDI {
14 using namespace GIDI;
15 #endif
16 
18 /*
19 ************************************************************
20 */
22 
23  MCGIDI_angularEnergy *angularEnergy;
24 
25  if( ( angularEnergy = (MCGIDI_angularEnergy *) smr_malloc2( smr, sizeof( MCGIDI_angularEnergy ), 0, "angularEnergy" ) ) == NULL ) return( NULL );
26  if( MCGIDI_angularEnergy_initialize( smr, angularEnergy ) ) angularEnergy = MCGIDI_angularEnergy_free( smr, angularEnergy );
27  return( angularEnergy );
28 }
29 /*
30 ************************************************************
31 */
33 
34  memset( angularEnergy, 0, sizeof( MCGIDI_angularEnergy ) );
35  return( 0 );
36 }
37 /*
38 ************************************************************
39 */
41 
42  MCGIDI_angularEnergy_release( smr, angularEnergy );
43  smr_freeMemory( (void **) &angularEnergy );
44  return( NULL );
45 }
46 /*
47 ************************************************************
48 */
50 
51  int i;
52 
53  for( i = 0; i < angularEnergy->pdfOfMuGivenE.numberOfWs; i++ ) {
54  MCGIDI_sampling_pdfsOfXGivenW_release( smr, &(angularEnergy->pdfOfEpGivenEAndMu[i]) );
55  }
56  smr_freeMemory( (void **) &(angularEnergy->pdfOfEpGivenEAndMu) );
57  MCGIDI_sampling_pdfsOfXGivenW_release( smr, &(angularEnergy->pdfOfMuGivenE) );
58  MCGIDI_angularEnergy_initialize( smr, angularEnergy );
59 
60  return( 0 );
61 }
62 /*
63 ************************************************************
64 */
66 
67  xDataTOM_element *angularEnergyElement, *pointwise = NULL;
68  char const *nativeData;
69 
70  if( ( angularEnergyElement = xDataTOME_getOneElementByName( smr, element, "angularEnergy", 1 ) ) == NULL ) goto err;
71 
72  if( ( nativeData = xDataTOM_getAttributesValueInElement( angularEnergyElement, "nativeData" ) ) == NULL ) goto err;
73  if( strcmp( nativeData, "pointwise" ) == 0 ) {
74  if( ( pointwise = xDataTOME_getOneElementByName( smr, angularEnergyElement, "pointwise", 1 ) ) == NULL ) goto err; }
75  else if( strcmp( nativeData, "linear" ) == 0 ) {
76  if( ( pointwise = xDataTOME_getOneElementByName( smr, angularEnergyElement, "linear", 1 ) ) == NULL ) goto err; }
77  else {
78  smr_setReportError2( smr, smr_unknownID, 1, "angularEnergy nativeData = '%s' not supported", nativeData );
79  goto err;
80  }
81  if( pointwise != NULL ) return( MCGIDI_angularEnergy_parsePointwiseFromTOM( smr, pointwise, distribution ) );
82 
83  return( 0 );
84 
85 err:
86  return( 1 );
87 }
88 /*
89 ************************************************************
90 */
92 
93  int iV, iW;
94  double y, norm, energyInFactor;
95  char const *energyUnit, *energyOutProbabilityUnits[2] = { "MeV", "1/MeV" };
96  MCGIDI_angularEnergy *angularEnergy = NULL;
97  ptwXY_interpolation interpolationXY, interpolationWY, interpolationVY;
98  xDataTOM_XYs *XYs;
99  xDataTOM_W_XYs *W_XYs;
100  xDataTOM_V_W_XYs *V_W_XYs;
101  MCGIDI_pdfsOfXGivenW *pdfOfMuGivenE, *pdfOfEpGivenEAndMu = NULL, *pdfOfEpGivenEAndMu2 = NULL;
102  ptwXYPoints *pdfXY1 = NULL, *pdfXY2 = NULL;
103  nfu_status status;
104 
105  if( MCGIDI_fromTOM_interpolation( smr, pointwise, 0, &interpolationVY ) ) goto err;
106  if( MCGIDI_fromTOM_interpolation( smr, pointwise, 1, &interpolationWY ) ) goto err;
107  if( MCGIDI_fromTOM_interpolation( smr, pointwise, 2, &interpolationXY ) ) goto err;
108  if( ( angularEnergy = MCGIDI_angularEnergy_new( smr ) ) == NULL ) goto err;
109 
110  if( ( angularEnergy->frame = MCGIDI_misc_getProductFrame( smr, pointwise ) ) == xDataTOM_frame_invalid ) goto err;
111 
112  pdfOfMuGivenE = &(angularEnergy->pdfOfMuGivenE);
113  pdfOfMuGivenE->interpolationWY = interpolationVY;
114  pdfOfMuGivenE->interpolationXY = interpolationWY;
115 
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_malloc2( smr, V_W_XYs->length * sizeof( double ), 1, "pdfOfMuGivenE->Ws" ) ) == NULL ) goto err;
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_pdfsOfXGivenW *) smr_malloc2( smr, V_W_XYs->length * sizeof( MCGIDI_pdfsOfXGivenW ), 1, "pdfOfEpGivenEAndMu" ) ) == NULL ) goto err;
120 
121  energyUnit = xDataTOM_subAxes_getUnit( smr, &(V_W_XYs->subAxes), 0 );
122  if( !smr_isOk( smr ) ) goto err;
123  energyInFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyUnit, "MeV" );
124  if( !smr_isOk( smr ) ) goto err;
125 
126  for( iV = 0; iV < V_W_XYs->length; iV++ ) {
127  W_XYs = &(V_W_XYs->W_XYs[iV]);
128  pdfOfEpGivenEAndMu2 = &(pdfOfEpGivenEAndMu[iV]);
129  pdfOfEpGivenEAndMu2->interpolationWY = interpolationWY;
130  pdfOfEpGivenEAndMu2->interpolationXY = interpolationXY;
131  if( ( pdfXY2 = ptwXY_new( interpolationWY, NULL, 2., 1e-6, W_XYs->length, 10, &status, 0 ) ) == NULL ) goto errA;
132  if( ( pdfOfEpGivenEAndMu2->Ws = (double *) smr_malloc2( smr, W_XYs->length * sizeof( double ), 1, "pdfOfEpGivenEAndMu2->Ws" ) ) == NULL ) goto err;
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++ ) {
135  XYs = &(W_XYs->XYs[iW]);
136  if( ( pdfXY1 = MCGIDI_misc_dataFromXYs2ptwXYPointsInUnitsOf( smr, XYs, interpolationXY, energyOutProbabilityUnits ) ) == NULL ) goto err;
137  y = ptwXY_integrateDomain( pdfXY1, &status );
138  if( ( status = ptwXY_setValueAtX( pdfXY2, XYs->value, y ) ) != nfu_Okay ) goto errA;
139  if( status != nfu_Okay ) goto errA;
140 
141  if( y == 0 ) {
142  if( ( status = ptwXY_add_double( pdfXY1, 0.5 ) ) != nfu_Okay ) goto errA;
143  }
144  pdfOfEpGivenEAndMu2->Ws[iW] = XYs->value;
145  if( MCGIDI_fromTOM_pdfOfX( smr, pdfXY1, &(pdfOfEpGivenEAndMu2->dist[iW]), &norm ) ) goto err;
146  pdfOfEpGivenEAndMu2->numberOfWs++;
147 
148  pdfXY1 = ptwXY_free( pdfXY1 );
149  }
150  pdfOfMuGivenE->Ws[iV] = energyInFactor * W_XYs->value;
151  if( MCGIDI_fromTOM_pdfOfX( smr, pdfXY2, &(pdfOfMuGivenE->dist[iV]), &norm ) ) goto err;
152  pdfOfMuGivenE->numberOfWs++;
153 
154  pdfXY2 = ptwXY_free( pdfXY2 );
155  }
156 
157  angularEnergy->pdfOfEpGivenEAndMu = pdfOfEpGivenEAndMu;
158  distribution->angularEnergy = angularEnergy;
160 
161  return( 0 );
162 
163 errA:
164  smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_integrateDomain err = %d: %s\n", status, nfu_statusMessage( status ) );
165 err:
166  if( pdfXY1 != NULL ) ptwXY_free( pdfXY1 );
167  if( pdfXY2 != NULL ) ptwXY_free( pdfXY2 );
168 /* Need to free pdfOfEpGivenEAndMu. */
169  if( angularEnergy != NULL ) MCGIDI_angularEnergy_free( smr, angularEnergy );
170  return( 1 );
171 }
172 /*
173 ************************************************************
174 */
176  MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
177 
178  int status = MCGIDI_sampling_doubleDistribution( smr, &(angularEnergy->pdfOfMuGivenE), angularEnergy->pdfOfEpGivenEAndMu, modes, decaySamplingInfo );
179 
180  decaySamplingInfo->frame = angularEnergy->frame;
181  return( status );
182 }
183 
184 #if defined __cplusplus
185 }
186 #endif
187