ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MCGIDI_KalbachMann.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MCGIDI_KalbachMann.cc
1 /*
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5 #include <string.h>
6 #include <cmath>
7 
8 #if defined __cplusplus
9 #include "G4Exp.hh"
10 #include "G4Log.hh"
11 #include "G4Pow.hh"
12 namespace GIDI {
13 using namespace GIDI;
14 #endif
15 const double C1 = 0.04, C2 = 1.8e-6/*, C3 = 6.7e-7*/;
16 /*
17 const double Et1 = 130., Et3 = 41.;
18 */
19 #if defined __cplusplus
20 }
21 #endif
22 
23 #include "MCGIDI_fromTOM.h"
24 #include "MCGIDI_misc.h"
25 #include "MCGIDI_private.h"
26 
27 #if defined __cplusplus
28 namespace GIDI {
29 using namespace GIDI;
30 #endif
31 
32 
33 static int MCGIDI_KalbachMann_parseFromTOM2( statusMessageReporting *smr, int dataPerEout, int index, xDataTOM_KalbachMannCoefficients *coefficientsXData,
34  double energyInFactor, double energyOutFactor, MCGIDI_KalbachMann *KalbachMann );
35 static double MCGIDI_KalbachMann_S_a_or_b( double Z_AB, double N_AB, double Z_C, double N_C, double I_ab );
36 /*
37 ************************************************************
38 */
39 MCGIDI_KalbachMann *MCGIDI_KalbachMann_new( statusMessageReporting *smr, ptwXY_interpolation interpolationWY,
40  ptwXY_interpolation interpolationXY ) {
41 
42  MCGIDI_KalbachMann *KalbachMann;
43 
44  if( ( KalbachMann = (MCGIDI_KalbachMann *) smr_malloc2( smr, sizeof( MCGIDI_KalbachMann ), 0, "KalbachMann" ) ) == NULL ) return( NULL );
45  if( MCGIDI_KalbachMann_initialize( smr, KalbachMann, interpolationWY, interpolationXY ) ) KalbachMann = MCGIDI_KalbachMann_free( smr, KalbachMann );
46  return( KalbachMann );
47 }
48 /*
49 ************************************************************
50 */
51 int MCGIDI_KalbachMann_initialize( statusMessageReporting * /*smr*/, MCGIDI_KalbachMann *KalbachMann, ptwXY_interpolation interpolationWY, ptwXY_interpolation interpolationXY ) {
52 
53  memset( KalbachMann, 0, sizeof( MCGIDI_KalbachMann ) );
54  KalbachMann->dists.interpolationWY = interpolationWY;
55  KalbachMann->dists.interpolationXY = interpolationXY;
56  return( 0 );
57 }
58 /*
59 ************************************************************
60 */
61 MCGIDI_KalbachMann *MCGIDI_KalbachMann_free( statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann ) {
62 
63  MCGIDI_KalbachMann_release( smr, KalbachMann );
64  smr_freeMemory( (void **) &KalbachMann );
65  return( NULL );
66 }
67 /*
68 ************************************************************
69 */
70 int MCGIDI_KalbachMann_release( statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann ) {
71 
72  int i;
73  MCGIDI_pdfsOfXGivenW *dists = &(KalbachMann->dists);
74 
75  for( i = 0; i < dists->numberOfWs; i++ ) {
76  smr_freeMemory( (void **) &(KalbachMann->ras[i].rs) );
77  smr_freeMemory( (void **) &(dists->dist[i].Xs) );
78  }
79  smr_freeMemory( (void **) &(KalbachMann->ras) );
80  smr_freeMemory( (void **) &(dists->Ws) );
81  smr_freeMemory( (void **) &(dists->dist) );
82 
83  MCGIDI_KalbachMann_initialize( smr, KalbachMann, ptwXY_interpolationLinLin, ptwXY_interpolationLinLin );
84  return( 0 );
85 }
86 /*
87 ************************************************************
88 */
89 int MCGIDI_KalbachMann_parseFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution ) {
90 
91  MCGIDI_KalbachMann *KalbachMann = NULL;
92  xDataTOM_element *KalbachMannElement;
93  int index, dataPerEout = 3;
94  double energyInFactor, energyOutFactor;
95  xDataTOM_xDataInfo *xDataInfo;
96  xDataTOM_KalbachMann *KalbachMannXData;
97  ptwXY_interpolation interpolationXY, interpolationWY;
98  char const *energyFromUnit, *energyToUnit = "MeV";
99 
100  MCGIDI_POP *productPOP = distribution->product->pop;
101  double productZ = productPOP->Z, productA = productPOP->A, productN = productA - productZ;
102  MCGIDI_target_heated *targetHeated = MCGIDI_product_getTargetHeated( smr, distribution->product );
103  MCGIDI_POP *projectilePOP = MCGIDI_target_heated_getPOPForProjectile( smr, targetHeated );
104  double projectileZ = projectilePOP->Z, projectileA = projectilePOP->A, projectileN = projectileA - projectileZ;
105  MCGIDI_POP *targetPOP = MCGIDI_target_heated_getPOPForTarget( smr, targetHeated );
106  double targetZ = targetPOP->Z, targetA = targetPOP->A, targetN = targetA - targetZ;
107  double Ia = 0., Ib = 0., Ma = -1, mb = -1;
108 
109  if( ( targetA == 0 ) && ( targetZ == 6 ) ) { /* Special case for C_000 evaluation. */
110  targetN = 6;
111  targetA = 12;
112  }
113  if( ( KalbachMannElement = xDataTOME_getOneElementByName( smr, element, "KalbachMann", 1 ) ) == NULL ) goto err;
114 
115  if( MCGIDI_fromTOM_interpolation( smr, KalbachMannElement, 0, &interpolationWY ) ) goto err;
116  if( MCGIDI_fromTOM_interpolation( smr, KalbachMannElement, 1, &interpolationXY ) ) goto err;
117 
118  xDataInfo = &(KalbachMannElement->xDataInfo);
119  KalbachMannXData = (xDataTOM_KalbachMann *) xDataInfo->data;
120  if( KalbachMannXData->type == xDataTOM_KalbachMannType_fra ) dataPerEout = 4;
121 
122  energyFromUnit = xDataTOM_axes_getUnit( smr, &(xDataInfo->axes), 0 );
123  if( !smr_isOk( smr ) ) goto err;
124  energyInFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyFromUnit, energyToUnit );
125  if( !smr_isOk( smr ) ) goto err;
126 
127  energyFromUnit = xDataTOM_axes_getUnit( smr, &(xDataInfo->axes), 1 );
128  if( !smr_isOk( smr ) ) goto err;
129  energyOutFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyFromUnit, energyToUnit );
130  if( !smr_isOk( smr ) ) goto err;
131 
132  if( ( KalbachMann = distribution->KalbachMann = MCGIDI_KalbachMann_new( smr, interpolationWY, interpolationXY ) ) == NULL ) goto err;
133 
134 /*
135  double productMass MCGIDI_product_getMass_MeV( smr, distribution->product ), residualMass;
136 */
137  KalbachMann->energyToMeVFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyToUnit, "MeV" );
138  KalbachMann->massFactor = (double) productZ + productN; /* This is not correct as masses are needed not Z and N. */
139  KalbachMann->massFactor /= projectileN + projectileZ + targetZ + targetN - productZ + productN;
140  KalbachMann->massFactor += 1.;
141 
142  if( projectileZ == 0 ) {
143  if( projectileN == 1 ) Ma = 1; }
144  else if( projectileZ == 1 ) {
145  if( projectileN == 1 ) {
146  Ma = 1; }
147  else if( projectileN == 2 ) {
148  Ia = 2.22;
149  Ma = 1; } }
150  else if( projectileZ == 2 ) {
151  if( projectileN == 2 ) {
152  Ia = 28.3;
153  Ma = 0;
154  }
155  }
156 
157  if( productZ == 0 ) {
158  if( productN == 1 ) mb = 0.5; }
159  else if( productZ == 1 ) {
160  if( productN == 1 ) {
161  mb = 1; }
162  else if( productN == 2 ) {
163  Ia = 2.22;
164  mb = 1; }
165  else if( productN == 3 ) {
166  Ib = 8.48;
167  mb = 1; } }
168  else if( productZ == 2 ) {
169  if( productN == 1 ) {
170  Ib = 7.72;
171  mb = 1; }
172  else if( productN == 2 ) {
173  Ib = 28.3;
174  mb = 2;
175  }
176  }
177 
178  KalbachMann->Ma = Ma;
179  KalbachMann->mb = mb;
180 
181  KalbachMann->Sa = MCGIDI_KalbachMann_S_a_or_b( targetZ, targetN, targetZ + projectileZ, targetN + projectileN, Ia );
182  KalbachMann->Sb = MCGIDI_KalbachMann_S_a_or_b( projectileZ + targetZ - productZ, projectileN + targetN - productN,
183  targetZ + projectileZ, targetN + projectileN, Ib );
184 
185  KalbachMann->dists.numberOfWs = 0;
186  if( ( KalbachMann->dists.Ws = (double *) smr_malloc2( smr, KalbachMannXData->numberOfEnergies * sizeof( double ), 0, "KalbachMann->dists->Ws" ) ) == NULL ) goto err;
187  if( ( KalbachMann->dists.dist = (MCGIDI_pdfOfX *) smr_malloc2( smr, KalbachMannXData->numberOfEnergies * sizeof( MCGIDI_pdfOfX ), 0, "KalbachMann->dists->dist" ) ) == NULL ) goto err;
188  if( ( KalbachMann->ras = (MCGIDI_KalbachMann_ras *) smr_malloc2( smr, KalbachMannXData->numberOfEnergies * sizeof( MCGIDI_KalbachMann_ras ), 0, "KalbachMann->ras" ) ) == NULL ) goto err;
189 
190  for( index = 0; index < KalbachMannXData->numberOfEnergies; index++ ) {
191  if( MCGIDI_KalbachMann_parseFromTOM2( smr, dataPerEout, index, &(KalbachMannXData->coefficients[index]),
192  energyInFactor, energyOutFactor, KalbachMann ) ) goto err;
193  }
194 
195  if( ( KalbachMann->frame = MCGIDI_misc_getProductFrame( smr, KalbachMannElement ) ) == xDataTOM_frame_invalid ) goto err;
196  distribution->type = MCGIDI_distributionType_KalbachMann_e;
197 
198  return( 0 );
199 
200 err:
201  if( KalbachMann != NULL ) MCGIDI_KalbachMann_free( smr, KalbachMann );
202  return( 1 );
203 }
204 /*
205 ************************************************************
206 */
207 static int MCGIDI_KalbachMann_parseFromTOM2( statusMessageReporting *smr, int dataPerEout, int index, xDataTOM_KalbachMannCoefficients *coefficientsXData,
208  double energyInFactor, double energyOutFactor, MCGIDI_KalbachMann *KalbachMann ) {
209 
210  int i, j, n = coefficientsXData->length / dataPerEout;
211  MCGIDI_pdfsOfXGivenW *dists = &(KalbachMann->dists);
212  MCGIDI_pdfOfX *dist = &(dists->dist[index]);
213  double norm, *p, *rs = NULL, *as_ = NULL, *Xs = NULL, *pdf, *cdf;
214  nfu_status status;
215  ptwXYPoints *pdfXY = NULL;
216  ptwXYPoint *point;
217  ptwXPoints *cdfX = NULL;
218  char const *ptwFunc = "";
219 
220  if( ( Xs = (double *) smr_malloc2( smr, 3 * n * sizeof( double ), 0, "Xs" ) ) == NULL ) goto err;
221  pdf = &(Xs[n]);
222  cdf = &(pdf[n]);
223 
224  if( ( rs = (double *) smr_malloc2( smr, ( dataPerEout - 2 ) * n * sizeof( double ), 0, "rs" ) ) == NULL ) goto err;
225  if( dataPerEout == 4 ) as_ = &(rs[n]);
226 
227  ptwFunc = "ptwXY_new";
228  if( ( pdfXY = ptwXY_new( KalbachMann->dists.interpolationXY, NULL, 2., 1e-3, n, 10, &status, 0 ) ) == NULL ) goto errXY;
229 
230  ptwFunc = "ptwXY_setXYPairAtIndex";
231  for( i = 0, p = coefficientsXData->coefficients; i < n; i++, p += dataPerEout ) {
232  if( ( status = ptwXY_setValueAtX( pdfXY, p[0], p[1] ) ) != nfu_Okay ) goto errXY;
233  rs[i] = p[2];
234  if( dataPerEout == 4 ) as_[i] = p[3];
235  }
236 
237  for( j = 0; j < n; j++ ) {
238  point = ptwXY_getPointAtIndex_Unsafely( pdfXY, j );
239  Xs[j] = energyOutFactor * point->x;
240  pdf[j] = point->y / energyOutFactor;
241  }
242 
243  ptwFunc = "ptwXY_runningIntegral";
244  if( ( cdfX = ptwXY_runningIntegral( pdfXY, &status ) ) == NULL ) goto errXY;
245  norm = ptwX_getPointAtIndex_Unsafely( cdfX, n - 1 );
246  if( std::fabs( 1. - norm ) > 0.99 ) {
247  smr_setReportError2( smr, smr_unknownID, 1, "bad norm = %e for angular.linear data", norm );
248  goto err;
249  }
250  for( j = 0; j < n; j++ ) cdf[j] = ptwX_getPointAtIndex_Unsafely( cdfX, j ) / norm;
251  for( j = 0; j < n; j++ ) pdf[j] /= norm;
252 
253  dists->numberOfWs++;
254  dists->Ws[index] = energyInFactor * coefficientsXData->value;
255  dist->numberOfXs = n;
256  dist->Xs = Xs;
257  dist->pdf = pdf;
258  dist->cdf = cdf;
259  KalbachMann->ras[index].rs = rs;
260  KalbachMann->ras[index].as = as_;
261 
262  pdfXY = ptwXY_free( pdfXY );
263  cdfX = ptwX_free( cdfX );
264  return( 0 );
265 
266 errXY:
267  smr_setReportError2( smr, smr_unknownID, 1, "%s error = %d: %s\n", ptwFunc, status, nfu_statusMessage( status ) );
268 
269 err:
270  if( Xs != NULL ) smr_freeMemory( (void **) &Xs);
271  if( rs != NULL ) smr_freeMemory( (void **) &rs);
272  if( pdfXY != NULL ) ptwXY_free( pdfXY );
273  if( cdfX != NULL ) cdfX = ptwX_free( cdfX );
274  return( 1 );
275 }
276 /*
277 ************************************************************
278 */
279 static double MCGIDI_KalbachMann_S_a_or_b( double Z_AB, double N_AB, double Z_C, double N_C, double I_ab ) {
280 
281  double A_AB = Z_AB + N_AB, A_C = Z_C + N_C;
282  double invA_AB_third = 1.0 / G4Pow::GetInstance()->A13( A_AB ), invA_C_third = 1.0 / G4Pow::GetInstance()->A13 ( A_C );
283  double NZA_AB = ( N_AB - Z_AB ) * ( N_AB - Z_AB ) / A_AB, NZA_C = ( N_C - Z_C ) * ( N_C - Z_C ) / A_C, S;
284 
285  S = 15.68 * ( A_C - A_AB ) - 28.07 * ( NZA_C - NZA_AB )
286  - 18.56 * ( A_C * invA_C_third - A_AB * invA_AB_third ) + 33.22 * ( NZA_C * invA_C_third - NZA_AB * invA_AB_third )
287  - 0.717 * ( Z_C * Z_C * invA_C_third - Z_AB * Z_AB * invA_AB_third ) + 1.211 * ( Z_C * Z_C / A_C - Z_AB * Z_AB / A_AB )
288  - I_ab;
289  return( S );
290 }
291 /*
292 ************************************************************
293 */
294 int MCGIDI_KalbachMann_sampleEp( statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann, MCGIDI_quantitiesLookupModes &modes,
295  MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
296 
297  double Epl, Epu, Ep, r, r2, rl, ru, a, a2, al, au, mu, randomEp = decaySamplingInfo->rng( decaySamplingInfo->rngState );
298  MCGIDI_pdfsOfXGivenW_sampled sampled;
299  MCGIDI_pdfsOfXGivenW *dists = &(KalbachMann->dists);
300  ptwXY_interpolation interpolationWY;
301 
302  sampled.smr = smr;
303  sampled.w = modes.getProjectileEnergy( );
304  MCGIDI_sampling_sampleX_from_pdfsOfXGivenW( dists, &sampled, randomEp );
305 
306  interpolationWY = sampled.interpolationWY;
307  if( sampled.iW < 0 ) {
308  interpolationWY = ptwXY_interpolationFlat;
309  if( sampled.iW == -2 ) { /* ???????????? This should probably report a warning. */
310  sampled.iW = 0; }
311  else if( sampled.iW == -1 ) {
312  sampled.iW = dists->numberOfWs - 1;
313  }
314  }
315 
316  Ep = sampled.x; /* Sampled Ep. */
317  if( sampled.interpolationXY == ptwXY_interpolationFlat ) { /* Now sample r. */
318  r = KalbachMann->ras[sampled.iW].rs[sampled.iX1]; }
319  else {
320  Epl = dists->dist[sampled.iW].Xs[sampled.iX1];
321  Epu = dists->dist[sampled.iW].Xs[sampled.iX1+1];
322  rl = KalbachMann->ras[sampled.iW].rs[sampled.iX1];
323  ru = KalbachMann->ras[sampled.iW].rs[sampled.iX1+1];
324  r = ( ru - rl ) / ( Epu - Epl ) * ( Ep - Epl ) + rl;
325  }
326  if( interpolationWY == ptwXY_interpolationLinLin ) {
327  if( sampled.interpolationXY == ptwXY_interpolationFlat ) {
328  r2 = KalbachMann->ras[sampled.iW+1].rs[sampled.iX2]; }
329  else {
330  Epl = dists->dist[sampled.iW+1].Xs[sampled.iX2];
331  Epu = dists->dist[sampled.iW+1].Xs[sampled.iX2+1];
332  rl = KalbachMann->ras[sampled.iW+1].rs[sampled.iX2];
333  ru = KalbachMann->ras[sampled.iW+1].rs[sampled.iX2+1];
334  r2 = ( ru - rl ) / ( Epu - Epl ) * ( Ep - Epl ) + rl;
335  }
336  r = sampled.frac * r + ( 1. - sampled.frac ) * r2;
337  }
338 
339  if( KalbachMann->ras[0].as == NULL ) { /* Now determine a. */
340  double X1, X3_2;
341  double eb = KalbachMann->massFactor * KalbachMann->energyToMeVFactor * Ep + KalbachMann->Sb;
342 
343  X1 = eb; /* Not valid for ea > Et1. */
344  X3_2 = eb * eb; /* Not valid for ea > Et3. */
345  a = X1 * ( C1 + C2 * X1 * X1 ) + C2 * KalbachMann->Ma * KalbachMann->mb * X3_2 * X3_2; }
346  else {
347  if( sampled.interpolationXY == ptwXY_interpolationFlat ) {
348  a = KalbachMann->ras[sampled.iW].as[sampled.iX1]; }
349  else {
350  Epl = dists->dist[sampled.iW].Xs[sampled.iX1];
351  Epu = dists->dist[sampled.iW].Xs[sampled.iX1+1];
352  al = KalbachMann->ras[sampled.iW].as[sampled.iX1];
353  au = KalbachMann->ras[sampled.iW].as[sampled.iX1+1];
354  a = ( au - al ) / ( Epu - Epl ) * ( Ep - Epl ) + al;
355  }
356  a2 = 0.;
357  if( interpolationWY == ptwXY_interpolationLinLin ) {
358  if( sampled.interpolationXY == ptwXY_interpolationFlat ) {
359  a2 = KalbachMann->ras[sampled.iW+1].as[sampled.iX2]; }
360  else {
361  Epl = dists->dist[sampled.iW+1].Xs[sampled.iX2];
362  Epu = dists->dist[sampled.iW+1].Xs[sampled.iX2+1];
363  al = KalbachMann->ras[sampled.iW+1].as[sampled.iX2];
364  au = KalbachMann->ras[sampled.iW+1].as[sampled.iX2+1];
365  a2 = ( au - al ) / ( Epu - Epl ) * ( Ep - Epl ) + al;
366  }
367  }
368  a = sampled.frac * a + ( 1. - sampled.frac ) * a2;
369  }
370 
371  /* In the following: Cosh[ a mu ] + r Sinh[ a mu ] = ( 1 - r ) Cosh[ a mu ] + r ( Cosh[ a mu ] + Sinh[ a mu ] ). */
372  if( decaySamplingInfo->rng( decaySamplingInfo->rngState ) >= r ) { /* Sample the '( 1 - r ) Cosh[ a mu ]' term. */
373  double T = ( 2. * decaySamplingInfo->rng( decaySamplingInfo->rngState ) - 1. ) * std::sinh( a );
374 
375  mu = G4Log( T + std::sqrt( T * T + 1. ) ) / a; }
376  else { /* Sample the 'r ( Cosh[ a mu ] + Sinh[ a mu ]' term. */
377  double rng1 = decaySamplingInfo->rng( decaySamplingInfo->rngState ), exp_a = G4Exp( a );
378 
379  mu = G4Log( rng1 * exp_a + ( 1. - rng1 ) / exp_a ) / a;
380  }
381  if( mu < -1 ) {
382  mu = -1;}
383  else if( mu > 1 ) {
384  mu = 1;
385  }
386 
387  decaySamplingInfo->frame = KalbachMann->frame;
388  decaySamplingInfo->Ep = Ep;
389  decaySamplingInfo->mu = mu;
390  return( !smr_isOk( smr ) );
391 }
392 
393 #if defined __cplusplus
394 }
395 #endif