Geant4 Cross Reference |
1 /* 1 2 # <<BEGIN-copyright>> 3 # <<END-copyright>> 4 */ 5 6 #include <stdlib.h> 7 #include <cmath> 8 #include <float.h> 9 10 #include "ptwXY.h" 11 12 #if defined __cplusplus 13 #include <cmath> 14 #include "G4Exp.hh" 15 #include "G4Log.hh" 16 namespace GIDI { 17 using namespace GIDI; 18 #endif 19 20 static nfu_status ptwXY_createGaussianCentered 21 /* 22 ********************************************** 23 */ 24 ptwXPoints *ptwXY_getXArray( ptwXYPoints *ptwX 25 26 int64_t i, n; 27 ptwXPoints *xArray; 28 29 if( ( *status = ptwXY->status ) != nfu_Oka 30 n = ptwXY->length; 31 32 if( ( *status = ptwXY_simpleCoalescePoints 33 if( ( xArray = ptwX_new( n, status ) ) == 34 for( i = 0; i < n; i++ ) xArray->points[i] 35 xArray->length = n; 36 37 return( xArray ); 38 } 39 /* 40 ********************************************** 41 */ 42 nfu_status ptwXY_dullEdges( ptwXYPoints *ptwXY 43 44 #define minEps 5e-16 45 46 nfu_status status; 47 double xm, xp, dx, y, x1, y1, x2, y2, sign 48 ptwXYPoint *p; 49 50 /* This routine can only be used for linear in 51 This needs to be fixed and documented. */ 52 if( ( status = ptwXY->status ) != nfu_Okay 53 if( ptwXY->interpolation == ptwXY_interpol 54 if( ptwXY->interpolation == ptwXY_interpol 55 56 if( ptwXY->length < 2 ) return( nfu_Okay ) 57 58 if( lowerEps != 0. ) { 59 if( std::fabs( lowerEps ) < minEps ) { 60 sign = 1; 61 if( lowerEps < 0. ) sign = -1; 62 lowerEps = sign * minEps; 63 } 64 65 p = ptwXY_getPointAtIndex_Unsafely( pt 66 x1 = p->x; 67 y1 = p->y; 68 p = ptwXY_getPointAtIndex_Unsafely( pt 69 x2 = p->x; 70 y2 = p->y; 71 72 if( y1 != 0. ) { 73 dx = std::fabs( x1 * lowerEps ); 74 if( x1 == 0 ) dx = std::fabs( lowe 75 xm = x1 - dx; 76 xp = x1 + dx; 77 if( ( xp + dx ) < x2 ) { 78 if( ( status = ptwXY_getValueA 79 if( ( status = ptwXY_setValueA 80 else { 81 xp = x2; 82 y = y2; 83 } 84 if( lowerEps > 0 ) { 85 if( ( status = ptwXY_setValueA 86 else { 87 if( ( xm < 0. ) && ( x1 >= 0. 88 if( ( status = ptwXY_setVa 89 else { 90 if( ( status = ptwXY_setVa 91 if( ( status = ptwXY_inter 92 if( ( status = ptwXY_setVa 93 } 94 } 95 } 96 } 97 98 if( upperEps != 0. ) { 99 if( std::fabs( upperEps ) < minEps ) { 100 sign = 1; 101 if( upperEps < 0. ) sign = -1; 102 upperEps = sign * minEps; 103 } 104 105 p = ptwXY_getPointAtIndex_Unsafely( pt 106 x1 = p->x; 107 y1 = p->y; 108 p = ptwXY_getPointAtIndex_Unsafely( pt 109 x2 = p->x; 110 y2 = p->y; 111 112 if( y2 != 0. ) { 113 dx = std::fabs( x2 * upperEps ); 114 if( x2 == 0 ) dx = std::fabs( uppe 115 xm = x2 - dx; 116 xp = x2 + dx; 117 if( ( xm - dx ) > x1 ) { 118 if( ( status = ptwXY_getValueA 119 if( ( status = ptwXY_setValueA 120 else { 121 xm = x1; 122 y = y1; 123 } 124 if( upperEps < 0 ) { 125 if( ( status = ptwXY_setValueA 126 else { 127 if( ( status = ptwXY_setValueA 128 if( ( status = ptwXY_interpola 129 if( ( status = ptwXY_setValueA 130 } 131 } 132 } 133 134 return( ptwXY->status ); 135 136 #undef minEps 137 } 138 /* 139 ********************************************** 140 */ 141 nfu_status ptwXY_mergeClosePoints( ptwXYPoints 142 143 int64_t i, i1, j, k, n = ptwXY->length; 144 double x, y; 145 ptwXYPoint *p1, *p2; 146 147 if( n < 2 ) return( ptwXY->status ); 148 if( epsilon < 4 * DBL_EPSILON ) epsilon = 149 if( ptwXY_simpleCoalescePoints( ptwXY ) != 150 151 p2 = ptwXY->points; 152 x = p2->x; 153 for( i1 = 1, p2++; i1 < ( n - 1 ); i1++, p 154 if( ( p2->x - x ) > 0.5 * epsilon * ( 155 } 156 if( i1 != 1 ) { 157 for( i = i1, p1 = &(ptwXY->points[1]); 158 n = ptwXY->length = ptwXY->length - i1 159 } 160 161 p1 = &(ptwXY->points[n-1]); 162 x = p1->x; 163 for( i1 = n - 2, p1--; i1 > 0; i1--, p1-- 164 if( x - p1->x > 0.5 * epsilon * ( std: 165 } 166 if( i1 != ( n - 2 ) ) { 167 ptwXY->points[i1 + 1] = ptwXY->points[ 168 n = ptwXY->length = i1 + 2; 169 } 170 171 for( i = 1; i < n - 1; i++ ) { 172 p1 = &(ptwXY->points[i]); 173 x = p1->x; 174 y = p1->y; 175 for( j = i + 1, p2 = &(ptwXY->points[i 176 if( ( p2->x - p1->x ) > 0.5 * epsi 177 x += p2->x; 178 y += p2->y; 179 } 180 if( ( k = ( j - i ) ) > 1 ) { 181 p1->x = x / k; 182 p1->y = y / k; 183 for( p1 = &(ptwXY->points[i+1]); j 184 n -= ( k - 1 ); 185 } 186 } 187 ptwXY->length = n; 188 189 return( ptwXY->status ); 190 } 191 /* 192 ********************************************** 193 */ 194 ptwXYPoints *ptwXY_intersectionWith_ptwX( ptwX 195 196 int64_t i, i1, i2, lengthX = ptwX_length( 197 double x, y, xMin, xMax; 198 ptwXYPoints *n = NULL; 199 200 if( ( *status = ptwXY->status ) != nfu_Oka 201 if( ( *status = ptwX->status ) != nfu_Okay 202 if( ( *status = ptwXY_simpleCoalescePoints 203 *status = nfu_otherInterpolation; 204 if( ptwXY->interpolation == ptwXY_interpol 205 if( ( n = ptwXY_clone( ptwXY, status ) ) = 206 if( ptwXY->length == 0 ) return( n ); 207 xMin = ptwXY->points[0].x; 208 xMax = ptwXY->points[ptwXY->length - 1].x; 209 210 if( ( xMin >= ptwX->points[lengthX-1] ) || 211 n->length = 0; 212 return( n ); 213 } 214 215 for( i = 0; i < lengthX; i++ ) { /* 216 x = ptwX->points[i]; 217 if( x <= xMin ) continue; 218 if( x >= xMax ) break; 219 if( ( *status = ptwXY_getValueAtX( ptw 220 if( ( *status = ptwXY_setValueAtX( n, 221 } 222 if( ( *status = ptwXY_simpleCoalescePoints 223 224 i1 = 0; 225 i2 = n->length - 1; 226 if( lengthX > 0 ) { 227 x = ptwX->points[0]; 228 if( x > n->points[i1].x ) { 229 for( ; i1 < n->length; i1++ ) { 230 if( n->points[i1].x == x ) bre 231 } 232 } 233 234 x = ptwX->points[lengthX - 1]; 235 if( x < n->points[i2].x ) { 236 for( ; i2 > i1; i2-- ) { 237 if( n->points[i2].x == x ) bre 238 } 239 } 240 } 241 i2++; 242 243 if( i1 != 0 ) { 244 for( i = i1; i < i2; i++ ) n->points[i 245 } 246 n->length = i2 - i1; 247 248 return( n ); 249 250 Err: 251 ptwXY_free( n ); 252 return( NULL ); 253 } 254 /* 255 ********************************************** 256 */ 257 nfu_status ptwXY_areDomainsMutual( ptwXYPoints 258 259 nfu_status status; 260 int64_t n1 = ptwXY1->length, n2 = ptwXY2-> 261 ptwXYPoint *xy1, *xy2; 262 263 if( ( status = ptwXY1->status ) != nfu_Oka 264 if( ( status = ptwXY2->status ) != nfu_Oka 265 if( n1 == 0 ) return( nfu_empty ); 266 if( n2 == 0 ) return( nfu_empty ); 267 if( n1 < 2 ) { 268 status = nfu_tooFewPoints; } 269 else if( n2 < 2 ) { 270 status = nfu_tooFewPoints; } 271 else { 272 xy1 = ptwXY_getPointAtIndex_Unsafely( 273 xy2 = ptwXY_getPointAtIndex_Unsafely( 274 if( xy1->x < xy2->x ) { 275 if( xy2->y != 0. ) status = nfu_do 276 else if( xy1->x > xy2->x ) { 277 if( xy1->y != 0. ) status = nfu_do 278 } 279 280 if( status == nfu_Okay ) { 281 xy1 = ptwXY_getPointAtIndex_Unsafe 282 xy2 = ptwXY_getPointAtIndex_Unsafe 283 if( xy1->x < xy2->x ) { 284 if( xy1->y != 0. ) status = nf 285 else if( xy1->x > xy2->x ) { 286 if( xy2->y != 0. ) status = nf 287 } 288 } 289 } 290 return( status ); 291 } 292 /* 293 ********************************************** 294 */ 295 nfu_status ptwXY_tweakDomainsToMutualify( ptwX 296 297 nfu_status status; 298 int64_t n1 = ptwXY1->length, n2 = ptwXY2-> 299 double sum, diff; 300 ptwXYPoint *xy1, *xy2; 301 302 epsilon = std::fabs( epsilon ) + std::fabs 303 304 if( ( status = ptwXY1->status ) != nfu_Oka 305 if( ( status = ptwXY2->status ) != nfu_Oka 306 if( n1 == 0 ) return( nfu_empty ); 307 if( n2 == 0 ) return( nfu_empty ); 308 if( n1 < 2 ) { 309 status = nfu_tooFewPoints; } 310 else if( n2 < 2 ) { 311 status = nfu_tooFewPoints; } 312 else { 313 xy1 = ptwXY_getPointAtIndex_Unsafely( 314 xy2 = ptwXY_getPointAtIndex_Unsafely( 315 if( xy1->x < xy2->x ) { 316 if( xy2->y != 0. ) { 317 sum = std::fabs( xy1->x ) + st 318 diff = std::fabs( xy2->x - xy1 319 if( diff > epsilon * sum ) { 320 status = nfu_domainsNotMut 321 else { 322 xy1->x = xy2->x; 323 } 324 } } 325 else if( xy1->x > xy2->x ) { 326 if( xy1->y != 0. ) { 327 sum = std::fabs( xy1->x ) + st 328 diff = std::fabs( xy2->x - xy1 329 if( diff > epsilon * sum ) { 330 status = nfu_domainsNotMut 331 else { 332 xy2->x = xy1->x; 333 } 334 } 335 } 336 337 if( status == nfu_Okay ) { 338 xy1 = ptwXY_getPointAtIndex_Unsafe 339 xy2 = ptwXY_getPointAtIndex_Unsafe 340 if( xy1->x < xy2->x ) { 341 if( xy1->y != 0. ) { 342 sum = std::fabs( xy1->x ) 343 diff = std::fabs( xy2->x - 344 if( diff > epsilon * sum ) 345 status = nfu_domainsNo 346 else { 347 xy2->x = xy1->x; 348 } 349 } } 350 else if( xy1->x > xy2->x ) { 351 if( xy2->y != 0. ) { 352 sum = std::fabs( xy1->x ) 353 diff = std::fabs( xy2->x - 354 if( diff > epsilon * sum ) 355 status = nfu_domainsNo 356 else { 357 xy1->x = xy2->x; 358 } 359 } 360 } 361 } 362 } 363 return( status ); 364 } 365 /* 366 ********************************************** 367 */ 368 nfu_status ptwXY_mutualifyDomains( ptwXYPoints 369 ptwX 370 371 nfu_status status; 372 int64_t n1 = ptwXY1->length, n2 = ptwXY2-> 373 ptwXYPoint *xy1, *xy2; 374 375 switch( status = ptwXY_areDomainsMutual( p 376 case nfu_Okay : 377 case nfu_empty : 378 return( nfu_Okay ); 379 case nfu_domainsNotMutual : 380 break; 381 default : 382 return( status ); 383 } 384 385 if( ptwXY1->interpolation == ptwXY_interpo 386 if( ptwXY2->interpolation == ptwXY_interpo 387 if( ptwXY1->interpolation == ptwXY_interpo 388 if( ptwXY2->interpolation == ptwXY_interpo 389 390 xy1 = ptwXY_getPointAtIndex_Unsafely( ptwX 391 xy2 = ptwXY_getPointAtIndex_Unsafely( ptwX 392 if( xy1->x < xy2->x ) { 393 lowerEps1 = 0.; 394 if( xy2->y == 0. ) lowerEps2 = 0.; } 395 else if( xy1->x > xy2->x ) { 396 lowerEps2 = 0.; 397 if( xy1->y == 0. ) lowerEps1 = 0.; } 398 else { 399 lowerEps1 = lowerEps2 = 0.; 400 } 401 402 xy1 = ptwXY_getPointAtIndex_Unsafely( ptwX 403 xy2 = ptwXY_getPointAtIndex_Unsafely( ptwX 404 if( xy1->x < xy2->x ) { 405 upperEps2 = 0.; 406 if( xy1->y == 0. ) upperEps1 = 0.; } 407 else if( xy1->x > xy2->x ) { 408 upperEps1 = 0.; 409 if( xy2->y == 0. ) upperEps2 = 0.; } 410 else { 411 upperEps1 = upperEps2 = 0.; 412 } 413 414 if( ( lowerEps1 != 0. ) || ( upperEps1 != 415 if( ( status = ptwXY_dullEdges( ptwXY1 416 if( ( lowerEps2 != 0. ) || ( upperEps2 != 417 if( ( status = ptwXY_dullEdges( ptwXY2 418 419 return( status ); 420 } 421 /* 422 ********************************************** 423 */ 424 nfu_status ptwXY_copyToC_XY( ptwXYPoints *ptwX 425 426 int64_t i; 427 double *d = xys; 428 nfu_status status; 429 ptwXYPoint *pointFrom; 430 431 if( ptwXY->status != nfu_Okay ) return( pt 432 if( ( status = ptwXY_simpleCoalescePoints( 433 434 if( index1 < 0 ) index1 = 0; 435 if( index2 > ptwXY->length ) index2 = ptwX 436 if( index2 < index1 ) index2 = index1; 437 *numberOfPoints = index2 - index1; 438 if( allocatedSize < ( index2 - index1 ) ) 439 440 for( i = index1, pointFrom = ptwXY->points 441 *(d++) = pointFrom->x; 442 *(d++) = pointFrom->y; 443 } 444 445 return( nfu_Okay ); 446 } 447 /* 448 ********************************************** 449 */ 450 nfu_status ptwXY_valueTo_ptwXAndY( ptwXYPoints 451 452 int64_t i1, length = ptwXY_length( ptwXY ) 453 double *xps, *yps; 454 ptwXYPoint *pointFrom; 455 nfu_status status; 456 457 if( ptwXY->status != nfu_Okay ) return( pt 458 if( ( status = ptwXY_simpleCoalescePoints( 459 460 if( ( *xs = (double *) malloc( length * si 461 if( ( *ys = (double *) malloc( length * si 462 free( *xs ); 463 *xs = NULL; 464 return( nfu_mallocError ); 465 } 466 467 for( i1 = 0, pointFrom = ptwXY->points, xp 468 *xps = pointFrom->x; 469 *yps = pointFrom->y; 470 } 471 472 return( nfu_Okay ); 473 } 474 /* 475 ********************************************** 476 */ 477 ptwXYPoints *ptwXY_valueTo_ptwXY( double x1, d 478 479 ptwXYPoints *n; 480 481 *status = nfu_XNotAscending; 482 if( x1 >= x2 ) return( NULL ); 483 *status = nfu_Okay; 484 if( ( n = ptwXY_new( ptwXY_interpolationLi 485 ptwXY_setValueAtX( n, x1, y ); 486 ptwXY_setValueAtX( n, x2, y ); 487 return( n ); 488 } 489 /* 490 ********************************************** 491 */ 492 ptwXYPoints *ptwXY_createGaussianCenteredSigma 493 494 int64_t i, n; 495 ptwXYPoint *pm, *pp; 496 double x1, y1, x2, y2, accuracy2, yMin = 1 497 ptwXYPoints *gaussian; 498 499 if( accuracy < 1e-5 ) accuracy = 1e-5; 500 if( accuracy > 1e-1 ) accuracy = 1e-1; 501 if( ( gaussian = ptwXY_new( ptwXY_interpol 502 accuracy2 = accuracy = gaussian->accuracy; 503 if( accuracy2 > 5e-3 ) accuracy2 = 5e-3; 504 505 x1 = -std::sqrt( -2. * G4Log( yMin ) ); 506 y1 = yMin; 507 x2 = -5.2; 508 y2 = G4Exp( -0.5 * x2 * x2 ); 509 if( ( *status = ptwXY_setValueAtX( gaussia 510 gaussian->accuracy = 20 * accuracy2; 511 if( ( *status = ptwXY_createGaussianCenter 512 x1 = x2; 513 y1 = y2; 514 x2 = -4.; 515 y2 = G4Exp( -0.5 * x2 * x2 ); 516 gaussian->accuracy = 5 * accuracy2; 517 if( ( *status = ptwXY_createGaussianCenter 518 x1 = x2; 519 y1 = y2; 520 x2 = -1; 521 y2 = G4Exp( -0.5 * x2 * x2 ); 522 gaussian->accuracy = accuracy; 523 if( ( *status = ptwXY_createGaussianCenter 524 x1 = x2; 525 y1 = y2; 526 x2 = 0; 527 y2 = G4Exp( -0.5 * x2 * x2 ); 528 if( ( *status = ptwXY_createGaussianCenter 529 530 n = gaussian->length; 531 if( ( *status = ptwXY_coalescePoints( gaus 532 if( ( *status = ptwXY_setValueAtX( gaussia 533 pp = &(gaussian->points[gaussian->length]) 534 for( i = 0, pm = pp - 2; i < n; i++, pp++, 535 *pp = *pm; 536 pp->x *= -1; 537 } 538 gaussian->length = 2 * n + 1; 539 540 return( gaussian ); 541 542 Err: 543 ptwXY_free( gaussian ); 544 return( NULL ); 545 } 546 /* 547 ********************************************** 548 */ 549 static nfu_status ptwXY_createGaussianCentered 550 551 nfu_status status = nfu_Okay; 552 int morePoints = 0; 553 double x = 0.5 * ( x1 + x2 ); 554 double y = G4Exp( -x * x / 2 ), yMin = ( y 555 556 if( std::fabs( y - yMin ) > y * ptwXY->acc 557 if( morePoints && ( status = ptwXY_createG 558 if( ( status = ptwXY_setValueAtX( ptwXY, x 559 if( morePoints && ( status = ptwXY_createG 560 if( addX1Point ) status = ptwXY_setValueAt 561 return( status ); 562 } 563 /* 564 ********************************************** 565 */ 566 ptwXYPoints *ptwXY_createGaussian( double accu 567 double /*dullEps*/, nfu_status *status 568 569 int64_t i; 570 ptwXYPoints *gaussian, *sliced; 571 ptwXYPoint *point; 572 573 if( ( gaussian = ptwXY_createGaussianCente 574 for( i = 0, point = gaussian->points; i < 575 point->x = point->x * sigma + xCenter; 576 point->y *= amplitude; 577 } 578 if( ( gaussian->points[0].x < xMin ) || ( 579 if( ( sliced = ptwXY_xSlice( gaussian, 580 ptwXY_free( gaussian ); 581 gaussian = sliced; 582 } 583 584 return( gaussian ); 585 586 Err: 587 ptwXY_free( gaussian ); 588 return( NULL ); 589 } 590 591 #if defined __cplusplus 592 } 593 #endif 594