Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MCGIDI_KalbachMann.cc File Reference
#include <string.h>
#include <cmath>
#include "MCGIDI_fromTOM.h"
#include "MCGIDI_misc.h"
#include "MCGIDI_private.h"
Include dependency graph for MCGIDI_KalbachMann.cc:

Go to the source code of this file.

Functions

static int MCGIDI_KalbachMann_parseFromTOM2 (statusMessageReporting *smr, int dataPerEout, int index, xDataTOM_KalbachMannCoefficients *coefficientsXData, double energyInFactor, double energyOutFactor, MCGIDI_KalbachMann *KalbachMann)
 
static double MCGIDI_KalbachMann_S_a_or_b (double Z_AB, double N_AB, double Z_C, double N_C, double I_ab)
 
MCGIDI_KalbachMannMCGIDI_KalbachMann_new (statusMessageReporting *smr, ptwXY_interpolation interpolationWY, ptwXY_interpolation interpolationXY)
 
int MCGIDI_KalbachMann_initialize (statusMessageReporting *, MCGIDI_KalbachMann *KalbachMann, ptwXY_interpolation interpolationWY, ptwXY_interpolation interpolationXY)
 
MCGIDI_KalbachMannMCGIDI_KalbachMann_free (statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann)
 
int MCGIDI_KalbachMann_release (statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann)
 
int MCGIDI_KalbachMann_parseFromTOM (statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution)
 
int MCGIDI_KalbachMann_sampleEp (statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)
 

Variables

const double C1 = 0.04
 
const double C2 = 1.8e-6
 

Function Documentation

MCGIDI_KalbachMann* MCGIDI_KalbachMann_free ( statusMessageReporting smr,
MCGIDI_KalbachMann KalbachMann 
)

Definition at line 61 of file MCGIDI_KalbachMann.cc.

61  {
62 
63  MCGIDI_KalbachMann_release( smr, KalbachMann );
64  smr_freeMemory( (void **) &KalbachMann );
65  return( NULL );
66 }
int MCGIDI_KalbachMann_release(statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann)
void * smr_freeMemory(void **p)

Here is the call graph for this function:

Here is the caller graph for this function:

int MCGIDI_KalbachMann_initialize ( statusMessageReporting ,
MCGIDI_KalbachMann KalbachMann,
ptwXY_interpolation  interpolationWY,
ptwXY_interpolation  interpolationXY 
)

Definition at line 51 of file MCGIDI_KalbachMann.cc.

51  {
52 
53  memset( KalbachMann, 0, sizeof( MCGIDI_KalbachMann ) );
54  KalbachMann->dists.interpolationWY = interpolationWY;
55  KalbachMann->dists.interpolationXY = interpolationXY;
56  return( 0 );
57 }
MCGIDI_pdfsOfXGivenW dists
Definition: MCGIDI.h:376
ptwXY_interpolation interpolationXY
Definition: MCGIDI.h:306
ptwXY_interpolation interpolationWY
Definition: MCGIDI.h:306

Here is the caller graph for this function:

MCGIDI_KalbachMann* MCGIDI_KalbachMann_new ( statusMessageReporting smr,
ptwXY_interpolation  interpolationWY,
ptwXY_interpolation  interpolationXY 
)

Definition at line 39 of file MCGIDI_KalbachMann.cc.

40  {
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 }
MCGIDI_KalbachMann * MCGIDI_KalbachMann_free(statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann)
#define smr_malloc2(smr, size, zero, forItem)
int MCGIDI_KalbachMann_initialize(statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann, ptwXY_interpolation interpolationWY, ptwXY_interpolation interpolationXY)

Here is the call graph for this function:

Here is the caller graph for this function:

int MCGIDI_KalbachMann_parseFromTOM ( statusMessageReporting smr,
xDataTOM_element element,
MCGIDI_distribution distribution 
)

Definition at line 89 of file MCGIDI_KalbachMann.cc.

89  {
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;
197 
198  return( 0 );
199 
200 err:
201  if( KalbachMann != NULL ) MCGIDI_KalbachMann_free( smr, KalbachMann );
202  return( 1 );
203 }
MCGIDI_pdfsOfXGivenW dists
Definition: MCGIDI.h:376
char const * xDataTOM_axes_getUnit(statusMessageReporting *smr, xDataTOM_axes *axes, int index)
MCGIDI_POP * MCGIDI_target_heated_getPOPForTarget(statusMessageReporting *smr, MCGIDI_target_heated *target)
xDataTOM_KalbachMannCoefficients * coefficients
Definition: xDataTOM.h:141
double energyToMeVFactor
Definition: MCGIDI.h:375
MCGIDI_target_heated * MCGIDI_product_getTargetHeated(statusMessageReporting *smr, MCGIDI_product *product)
int MCGIDI_fromTOM_interpolation(statusMessageReporting *smr, xDataTOM_element *element, int index, enum ptwXY_interpolation_e *interpolation)
MCGIDI_KalbachMann * MCGIDI_KalbachMann_free(statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann)
MCGIDI_KalbachMann * MCGIDI_KalbachMann_new(statusMessageReporting *smr, ptwXY_interpolation interpolationWY, ptwXY_interpolation interpolationXY)
xDataTOM_element * xDataTOME_getOneElementByName(statusMessageReporting *smr, xDataTOM_element *element, char const *name, int required)
Definition: xDataTOM.cc:246
#define smr_malloc2(smr, size, zero, forItem)
int smr_isOk(statusMessageReporting *smr)
enum ptwXY_interpolation_e ptwXY_interpolation
MCGIDI_product * product
Definition: MCGIDI.h:381
MCGIDI_KalbachMann_ras * ras
Definition: MCGIDI.h:377
enum MCGIDI_distributionType type
Definition: MCGIDI.h:382
enum xDataTOM_frame frame
Definition: MCGIDI.h:374
MCGIDI_KalbachMann * KalbachMann
Definition: MCGIDI.h:387
xDataTOM_axes axes
Definition: xDataTOM.h:153
static double MCGIDI_KalbachMann_S_a_or_b(double Z_AB, double N_AB, double Z_C, double N_C, double I_ab)
enum xDataTOM_frame MCGIDI_misc_getProductFrame(statusMessageReporting *smr, xDataTOM_element *frameElement)
Definition: MCGIDI_misc.cc:315
double MCGIDI_misc_getUnitConversionFactor(statusMessageReporting *smr, char const *fromUnit, char const *toUnit)
Definition: MCGIDI_misc.cc:381
MCGIDI_POP * MCGIDI_target_heated_getPOPForProjectile(statusMessageReporting *smr, MCGIDI_target_heated *target)
xDataTOM_xDataInfo xDataInfo
Definition: xDataTOM.h:187
MCGIDI_POP * pop
Definition: MCGIDI.h:401
static int MCGIDI_KalbachMann_parseFromTOM2(statusMessageReporting *smr, int dataPerEout, int index, xDataTOM_KalbachMannCoefficients *coefficientsXData, double energyInFactor, double energyOutFactor, MCGIDI_KalbachMann *KalbachMann)
enum xDataTOM_KalbachMannType type
Definition: xDataTOM.h:138
MCGIDI_pdfOfX * dist
Definition: MCGIDI.h:308

Here is the call graph for this function:

Here is the caller graph for this function:

static int MCGIDI_KalbachMann_parseFromTOM2 ( statusMessageReporting smr,
int  dataPerEout,
int  index,
xDataTOM_KalbachMannCoefficients coefficientsXData,
double  energyInFactor,
double  energyOutFactor,
MCGIDI_KalbachMann KalbachMann 
)
static

Definition at line 207 of file MCGIDI_KalbachMann.cc.

208  {
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 }
MCGIDI_pdfsOfXGivenW dists
Definition: MCGIDI.h:376
ptwXY_interpolation interpolationXY
Definition: MCGIDI.h:306
double ptwX_getPointAtIndex_Unsafely(ptwXPoints *ptwX, int64_t index)
Definition: ptwX_core.cc:215
double * cdf
Definition: MCGIDI.h:301
const char * p
Definition: xmltok.h:285
ptwXYPoint * ptwXY_getPointAtIndex_Unsafely(ptwXYPoints *ptwXY, int64_t index)
Definition: ptwXY_core.cc:684
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
double * Xs
Definition: MCGIDI.h:299
#define smr_setReportError2(smr, libraryID, code, fmt,...)
ptwXPoints * ptwX_free(ptwXPoints *ptwX)
Definition: ptwX_core.cc:158
enum nfu_status_e nfu_status
#define smr_malloc2(smr, size, zero, forItem)
ptwXPoints * ptwXY_runningIntegral(ptwXYPoints *ptwXY, nfu_status *status)
const G4int n
#define smr_unknownID
double y
Definition: ptwXY.h:62
MCGIDI_KalbachMann_ras * ras
Definition: MCGIDI.h:377
ptwXYPoints * ptwXY_new(ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo const *interpolationOtherInfo, double biSectionMax, double accuracy, int64_t primarySize, int64_t secondarySize, nfu_status *status, int userFlag)
Definition: ptwXY_core.cc:29
int numberOfXs
Definition: MCGIDI.h:298
void * smr_freeMemory(void **p)
double x
Definition: ptwXY.h:62
double * pdf
Definition: MCGIDI.h:300
const char * nfu_statusMessage(nfu_status status)
Definition: nf_utilities.cc:76
MCGIDI_pdfOfX * dist
Definition: MCGIDI.h:308

Here is the call graph for this function:

Here is the caller graph for this function:

int MCGIDI_KalbachMann_release ( statusMessageReporting smr,
MCGIDI_KalbachMann KalbachMann 
)

Definition at line 70 of file MCGIDI_KalbachMann.cc.

70  {
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 
84  return( 0 );
85 }
MCGIDI_pdfsOfXGivenW dists
Definition: MCGIDI.h:376
double * Xs
Definition: MCGIDI.h:299
MCGIDI_KalbachMann_ras * ras
Definition: MCGIDI.h:377
void * smr_freeMemory(void **p)
int MCGIDI_KalbachMann_initialize(statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann, ptwXY_interpolation interpolationWY, ptwXY_interpolation interpolationXY)
MCGIDI_pdfOfX * dist
Definition: MCGIDI.h:308

Here is the call graph for this function:

Here is the caller graph for this function:

static double MCGIDI_KalbachMann_S_a_or_b ( double  Z_AB,
double  N_AB,
double  Z_C,
double  N_C,
double  I_ab 
)
static

Definition at line 279 of file MCGIDI_KalbachMann.cc.

279  {
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 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
double S(double temp)
G4double A13(G4double A) const
Definition: G4Pow.hh:132

Here is the call graph for this function:

Here is the caller graph for this function:

int MCGIDI_KalbachMann_sampleEp ( statusMessageReporting smr,
MCGIDI_KalbachMann KalbachMann,
MCGIDI_quantitiesLookupModes modes,
MCGIDI_decaySamplingInfo decaySamplingInfo 
)

Definition at line 294 of file MCGIDI_KalbachMann.cc.

295  {
296 
297  double Epl, Epu, Ep, r, r2, rl, ru, a, a2, al, au, mu, randomEp = decaySamplingInfo->rng( decaySamplingInfo->rngState );
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 }
const double C2
MCGIDI_pdfsOfXGivenW dists
Definition: MCGIDI.h:376
const double C1
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
double energyToMeVFactor
Definition: MCGIDI.h:375
ptwXY_interpolation interpolationWY
Definition: MCGIDI.h:313
double * Xs
Definition: MCGIDI.h:299
statusMessageReporting * smr
Definition: MCGIDI.h:312
int smr_isOk(statusMessageReporting *smr)
enum ptwXY_interpolation_e ptwXY_interpolation
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
MCGIDI_KalbachMann_ras * ras
Definition: MCGIDI.h:377
int MCGIDI_sampling_sampleX_from_pdfsOfXGivenW(MCGIDI_pdfsOfXGivenW *dists, MCGIDI_pdfsOfXGivenW_sampled *sampled, double r)
enum xDataTOM_frame frame
Definition: MCGIDI.h:374
double(* rng)(void *)
Definition: MCGIDI.h:258
enum xDataTOM_frame frame
Definition: MCGIDI.h:256
ptwXY_interpolation interpolationXY
Definition: MCGIDI.h:313
double getProjectileEnergy(void) const
Definition: MCGIDI.h:97
MCGIDI_pdfOfX * dist
Definition: MCGIDI.h:308

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

const double C1 = 0.04

Definition at line 15 of file MCGIDI_KalbachMann.cc.

const double C2 = 1.8e-6

Definition at line 15 of file MCGIDI_KalbachMann.cc.