ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MCGIDI_energyAngular.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MCGIDI_energyAngular.cc
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 
17 /*
18 ************************************************************
19 */
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 */
32 
33  memset( energyAngular, 0, sizeof( MCGIDI_energyAngular ) );
34  return( 0 );
35 }
36 /*
37 ************************************************************
38 */
40 
41  MCGIDI_energyAngular_release( smr, energyAngular );
42  smr_freeMemory( (void **) &energyAngular );
43  return( NULL );
44 }
45 /*
46 ************************************************************
47 */
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 */
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 */
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  if( status != nfu_Okay ) goto errA;
146 
147  if( y == 0 ) {
148  if( ( status = ptwXY_add_double( pdfXY1, 0.5 ) ) != nfu_Okay ) goto errA;
149  }
150  pdfOfMuGivenEAndEp2->Ws[iW] = energyOutFactor * XYs->value;
151  if( MCGIDI_fromTOM_pdfOfX( smr, pdfXY1, &(pdfOfMuGivenEAndEp2->dist[iW]), &norm ) ) goto err;
152  pdfOfMuGivenEAndEp2->numberOfWs++;
153 
154  pdfXY1 = ptwXY_free( pdfXY1 );
155  }
156  pdfOfEpGivenE->Ws[iV] = energyInFactor * W_XYs->value;
157  if( MCGIDI_fromTOM_pdfOfX( smr, pdfXY2, &(pdfOfEpGivenE->dist[iV]), &norm ) ) goto err;
158  pdfOfEpGivenE->numberOfWs++;
159 
160  pdfXY2 = ptwXY_free( pdfXY2 );
161  }
162  energyAngular->pdfOfMuGivenEAndEp = pdfOfMuGivenEAndEp;
163  distribution->energyAngular = energyAngular;
165 
166  return( 0 );
167 
168 errA:
169  smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_integrateDomain err = %d: %s\n", status, nfu_statusMessage( status ) );
170 err:
171  if( pdfXY1 != NULL ) ptwXY_free( pdfXY1 );
172  if( pdfXY2 != NULL ) ptwXY_free( pdfXY2 );
173  if( energyAngular != NULL ) MCGIDI_energyAngular_free( smr, energyAngular );
174 /* ????????? Need to free pdfOfMuGivenEAndEp, now may be handled by MCGIDI_energyAngular_free. Need to check. */
175  return( 1 );
176 }
177 /*
178 ************************************************************
179 */
181  MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
182 
183  double Ep;
184  MCGIDI_energyAngular *energyAngular = distribution->energyAngular;
185 
186  MCGIDI_sampling_doubleDistribution( smr, &(energyAngular->pdfOfEpGivenE), energyAngular->pdfOfMuGivenEAndEp, modes, decaySamplingInfo );
187  Ep = decaySamplingInfo->mu;
188  decaySamplingInfo->mu = decaySamplingInfo->Ep;
189  decaySamplingInfo->Ep = Ep;
190  decaySamplingInfo->frame = energyAngular->frame;
191 
192  return( 0 );
193 }
194 
195 #if defined __cplusplus
196 }
197 #endif
198