Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // 28 // ------------------------------------------- 29 // first version created by P.Urban , 06/04/19 30 // modifications + "precise" functions added b 31 // modifications , TOF functions , 26/10/98, L 32 // cache mechanism in order to gain time, 11/0 33 // bug fixed , 12/04/99 , L.Urban 34 // 10.11.99: moved from RWT hash dictionary to 35 // 27.09.01 L.Urban , bug fixed (negative ener 36 // 26.10.01 all static functions moved from .i 37 // 15.01.03 Add interfaces required for "cut p 38 // 12.03.03 Add warnings to obsolete interface 39 // 10.04.03 Add call to G4LossTableManager is 40 // 41 // ------------------------------------------- 42 43 #include "G4EnergyLossTables.hh" 44 #include "G4SystemOfUnits.hh" 45 #include "G4MaterialCutsCouple.hh" 46 #include "G4RegionStore.hh" 47 #include "G4LossTableManager.hh" 48 49 50 //....oooOO0OOooo........oooOO0OOooo........oo 51 52 G4EnergyLossTablesHelper *G4EnergyLossTables:: 53 G4EnergyLossTablesHelper *G4EnergyLossTables:: 54 G4ParticleDefinition* G4EnergyLossTables::last 55 G4double G4EnergyLossTables::QQPositron = 1.0; 56 G4double G4EnergyLossTables::Chargesquare ; 57 G4int G4EnergyLossTables::oldIndex = -1 ; 58 G4double G4EnergyLossTables::rmin = 0. ; 59 G4double G4EnergyLossTables::rmax = 0. ; 60 G4double G4EnergyLossTables::Thigh = 0. ; 61 G4int G4EnergyLossTables::let_counter = 0; 62 G4int G4EnergyLossTables::let_max_num_warni 63 G4bool G4EnergyLossTables::first_loss = true 64 65 G4EnergyLossTables::helper_map *G4EnergyLossTa 66 67 //....oooOO0OOooo........oooOO0OOooo........oo 68 69 G4EnergyLossTablesHelper::G4EnergyLossTablesHe 70 const G4PhysicsTable* aDEDXTable, 71 const G4PhysicsTable* aRangeTable, 72 const G4PhysicsTable* anInverseRangeTable, 73 const G4PhysicsTable* aLabTimeTable, 74 const G4PhysicsTable* aProperTimeTable, 75 G4double aLowestKineticEnergy, 76 G4double aHighestKineticEnergy, 77 G4double aMassRatio, 78 G4int aNumberOfBins) 79 : 80 theDEDXTable(aDEDXTable), theRangeTable(aRan 81 theInverseRangeTable(anInverseRangeTable), 82 theLabTimeTable(aLabTimeTable), 83 theProperTimeTable(aProperTimeTable), 84 theLowestKineticEnergy(aLowestKineticEnergy) 85 theHighestKineticEnergy(aHighestKineticEnerg 86 theMassRatio(aMassRatio), 87 theNumberOfBins(aNumberOfBins) 88 { } 89 90 //....oooOO0OOooo........oooOO0OOooo........oo 91 92 G4EnergyLossTablesHelper::G4EnergyLossTablesHe 93 { 94 theLowestKineticEnergy = 0.0; 95 theHighestKineticEnergy= 0.0; 96 theMassRatio = 0.0; 97 theNumberOfBins = 0; 98 theDEDXTable = theRangeTable = theInverseRan 99 = theProperTimeTable = nullptr; 100 } 101 102 //....oooOO0OOooo........oooOO0OOooo........oo 103 104 void G4EnergyLossTables::Register( 105 const G4ParticleDefinition* p, 106 const G4PhysicsTable* tDEDX, 107 const G4PhysicsTable* tRange, 108 const G4PhysicsTable* tInverseRange, 109 const G4PhysicsTable* tLabTime, 110 const G4PhysicsTable* tProperTime, 111 G4double lowestKineticEnergy, 112 G4double highestKineticEnergy, 113 G4double massRatio, 114 G4int NumberOfBins) 115 { 116 if (!dict) dict = new G4EnergyLossTables::he 117 if (!null_loss) null_loss = new G4EnergyLoss 118 if (!t) t = new G4EnergyLossTablesHelper; 119 120 (*dict)[p]= G4EnergyLossTablesHelper(tDEDX, 121 tLabTime,tProperTime,lowes 122 highestKineticEnergy, massRatio,Number 123 124 *t = GetTables(p) ; // important for cach 125 lastParticle = (G4ParticleDefinition*) p ; 126 Chargesquare = (p->GetPDGCharge())*(p->GetPD 127 QQPositron ; 128 if (first_loss ) { 129 *null_loss = G4EnergyLossTablesHelper( 130 nullptr, nullptr, nullptr, nu 131 first_loss = false; 132 } 133 } 134 135 //....oooOO0OOooo........oooOO0OOooo........oo 136 137 const G4PhysicsTable* G4EnergyLossTables::GetD 138 const G4ParticleDefinition* p) 139 { 140 if (!dict) dict = new G4EnergyLossTables::he 141 helper_map::iterator it; 142 if((it=dict->find(p))==dict->end()) return n 143 return (*it).second.theDEDXTable; 144 } 145 146 //....oooOO0OOooo........oooOO0OOooo........oo 147 148 const G4PhysicsTable* G4EnergyLossTables::GetR 149 const G4ParticleDefinition* p) 150 { 151 if (!dict) dict = new G4EnergyLossTables::he 152 helper_map::iterator it; 153 if((it=dict->find(p))==dict->end()) return n 154 return (*it).second.theRangeTable; 155 } 156 157 //....oooOO0OOooo........oooOO0OOooo........oo 158 159 const G4PhysicsTable* G4EnergyLossTables::GetI 160 const G4ParticleDefinition* p) 161 { 162 if (!dict) dict = new G4EnergyLossTables::he 163 helper_map::iterator it; 164 if((it=dict->find(p))==dict->end()) return n 165 return (*it).second.theInverseRangeTable; 166 } 167 168 //....oooOO0OOooo........oooOO0OOooo........oo 169 170 const G4PhysicsTable* G4EnergyLossTables::GetL 171 const G4ParticleDefinition* p) 172 { 173 if (!dict) dict = new G4EnergyLossTables::he 174 helper_map::iterator it; 175 if((it=dict->find(p))==dict->end()) return n 176 return (*it).second.theLabTimeTable; 177 } 178 179 //....oooOO0OOooo........oooOO0OOooo........oo 180 181 const G4PhysicsTable* G4EnergyLossTables::GetP 182 const G4ParticleDefinition* p) 183 { 184 if (!dict) dict = new G4EnergyLossTables::he 185 helper_map::iterator it; 186 if((it=dict->find(p))==dict->end()) return n 187 return (*it).second.theProperTimeTable; 188 } 189 190 //....oooOO0OOooo........oooOO0OOooo........oo 191 192 G4EnergyLossTablesHelper G4EnergyLossTables::G 193 const G4ParticleDefinition* p) 194 { 195 if (!dict) dict = new G4EnergyLossTables::he 196 if (!null_loss) null_loss = new G4EnergyLoss 197 198 helper_map::iterator it; 199 if ((it=dict->find(p))==dict->end()) { 200 return *null_loss; 201 } 202 return (*it).second; 203 } 204 205 //....oooOO0OOooo........oooOO0OOooo........oo 206 207 G4double G4EnergyLossTables::GetDEDX( 208 const G4ParticleDefinition *aParticle, 209 G4double KineticEnergy, 210 const G4Material *aMaterial) 211 { 212 if (!t) t = new G4EnergyLossTablesHelper; 213 214 CPRWarning(); 215 if(aParticle != (const G4ParticleDefinition* 216 { 217 *t= GetTables(aParticle); 218 lastParticle = (G4ParticleDefinition*) aPa 219 Chargesquare = (aParticle->GetPDGCharge()) 220 (aParticle->GetPDGCharge()) 221 QQPositron ; 222 oldIndex = -1 ; 223 } 224 const G4PhysicsTable* dEdxTable= t->theDEDX 225 if (!dEdxTable) { 226 ParticleHaveNoLoss(aParticle,"dEdx"); 227 return 0.0; 228 } 229 230 G4int materialIndex = (G4int)aMaterial->GetI 231 G4double scaledKineticEnergy = KineticEnergy 232 G4double dEdx; 233 G4bool isOut; 234 235 if (scaledKineticEnergy<t->theLowestKineticE 236 237 dEdx =(*dEdxTable)(materialIndex)->GetVal 238 t->theLowestKineticEnergy,isOut) 239 *std::sqrt(scaledKineticEnergy/t->t 240 241 } else if (scaledKineticEnergy>t->theHighest 242 243 dEdx = (*dEdxTable)(materialIndex)->GetVa 244 t->theHighestKineticEnergy,isOut); 245 246 } else { 247 248 dEdx = (*dEdxTable)(materialIndex)->GetVal 249 scaledKineticEnergy,isOut); 250 251 } 252 253 return dEdx*Chargesquare; 254 } 255 256 //....oooOO0OOooo........oooOO0OOooo........oo 257 258 G4double G4EnergyLossTables::GetLabTime( 259 const G4ParticleDefinition *aParticle, 260 G4double KineticEnergy, 261 const G4Material *aMaterial) 262 { 263 if (!t) t = new G4EnergyLossTablesHelper; 264 265 CPRWarning(); 266 if(aParticle != (const G4ParticleDefinition* 267 { 268 *t= GetTables(aParticle); 269 lastParticle = (G4ParticleDefinition*) aPa 270 oldIndex = -1 ; 271 } 272 const G4PhysicsTable* labtimeTable= t->theLa 273 if (!labtimeTable) { 274 ParticleHaveNoLoss(aParticle,"LabTime"); 275 return 0.0; 276 } 277 278 const G4double parlowen=0.4 , ppar=0.5-parlo 279 G4int materialIndex = (G4int)aMaterial->GetI 280 G4double scaledKineticEnergy = KineticEnergy 281 G4double time; 282 G4bool isOut; 283 284 if (scaledKineticEnergy<t->theLowestKineticE 285 286 time = std::exp(ppar*std::log(scaledKinet 287 (*labtimeTable)(materialIndex)->Ge 288 t->theLowestKineticEnergy,isOut) 289 290 291 } else if (scaledKineticEnergy>t->theHighest 292 293 time = (*labtimeTable)(materialIndex)->Ge 294 t->theHighestKineticEnergy,isOut 295 296 } else { 297 298 time = (*labtimeTable)(materialIndex)->Get 299 scaledKineticEnergy,isOut); 300 301 } 302 303 return time/t->theMassRatio ; 304 } 305 306 //....oooOO0OOooo........oooOO0OOooo........oo 307 308 G4double G4EnergyLossTables::GetDeltaLabTime( 309 const G4ParticleDefinition *aParticle, 310 G4double KineticEnergyStart, 311 G4double KineticEnergyEnd, 312 const G4Material *aMaterial) 313 { 314 if (!t) t = new G4EnergyLossTablesHelper; 315 316 CPRWarning(); 317 if(aParticle != (const G4ParticleDefinition* 318 { 319 *t= GetTables(aParticle); 320 lastParticle = (G4ParticleDefinition*) aPa 321 oldIndex = -1 ; 322 } 323 const G4PhysicsTable* labtimeTable= t->theLa 324 if (!labtimeTable) { 325 ParticleHaveNoLoss(aParticle,"LabTime"); 326 return 0.0; 327 } 328 329 const G4double parlowen=0.4 , ppar=0.5-parlo 330 const G4double dToverT = 0.05 , facT = 1. -d 331 G4double timestart,timeend,deltatime,dTT; 332 G4bool isOut; 333 334 G4int materialIndex = (G4int)aMaterial->GetI 335 G4double scaledKineticEnergy = KineticEnergy 336 337 if (scaledKineticEnergy<t->theLowestKineticE 338 339 timestart = std::exp(ppar*std::log(scaled 340 (*labtimeTable)(materialIndex) 341 t->theLowestKineticEnergy,isOu 342 343 344 } else if (scaledKineticEnergy>t->theHighest 345 346 timestart = (*labtimeTable)(materialIndex 347 t->theHighestKineticEnergy,isO 348 349 } else { 350 351 timestart = (*labtimeTable)(materialIndex) 352 scaledKineticEnergy,isOut); 353 354 } 355 356 dTT = (KineticEnergyStart - KineticEnergyEnd 357 358 if( dTT < dToverT ) 359 scaledKineticEnergy = facT*KineticEnergySt 360 else 361 scaledKineticEnergy = KineticEnergyEnd*t-> 362 363 if (scaledKineticEnergy<t->theLowestKineticE 364 365 timeend = std::exp(ppar*std::log(scaledKi 366 (*labtimeTable)(materialIndex) 367 t->theLowestKineticEnergy,isOu 368 369 370 } else if (scaledKineticEnergy>t->theHighest 371 372 timeend = (*labtimeTable)(materialIndex)- 373 t->theHighestKineticEnergy,isO 374 375 } else { 376 377 timeend = (*labtimeTable)(materialIndex)-> 378 scaledKineticEnergy,isOut); 379 380 } 381 382 deltatime = timestart - timeend ; 383 384 if( dTT < dToverT ) 385 deltatime *= dTT/dToverT; 386 387 return deltatime/t->theMassRatio ; 388 } 389 390 //....oooOO0OOooo........oooOO0OOooo........oo 391 392 G4double G4EnergyLossTables::GetProperTime( 393 const G4ParticleDefinition *aParticle, 394 G4double KineticEnergy, 395 const G4Material *aMaterial) 396 { 397 if (!t) t = new G4EnergyLossTablesHelper; 398 399 CPRWarning(); 400 if(aParticle != (const G4ParticleDefinition* 401 { 402 *t= GetTables(aParticle); 403 lastParticle = (G4ParticleDefinition*) aPa 404 oldIndex = -1 ; 405 } 406 const G4PhysicsTable* propertimeTable= t->th 407 if (!propertimeTable) { 408 ParticleHaveNoLoss(aParticle,"ProperTime") 409 return 0.0; 410 } 411 412 const G4double parlowen=0.4 , ppar=0.5-parlo 413 G4int materialIndex = (G4int)aMaterial->GetI 414 G4double scaledKineticEnergy = KineticEnergy 415 G4double time; 416 G4bool isOut; 417 418 if (scaledKineticEnergy<t->theLowestKineticE 419 420 time = std::exp(ppar*std::log(scaledKinet 421 (*propertimeTable)(materialIndex)- 422 t->theLowestKineticEnergy,isOut) 423 424 425 } else if (scaledKineticEnergy>t->theHighest 426 427 time = (*propertimeTable)(materialIndex)- 428 t->theHighestKineticEnergy,isOut 429 430 } else { 431 432 time = (*propertimeTable)(materialIndex)-> 433 scaledKineticEnergy,isOut); 434 435 } 436 437 return time/t->theMassRatio ; 438 } 439 440 //....oooOO0OOooo........oooOO0OOooo........oo 441 442 G4double G4EnergyLossTables::GetDeltaProperTim 443 const G4ParticleDefinition *aParticle, 444 G4double KineticEnergyStart, 445 G4double KineticEnergyEnd, 446 const G4Material *aMaterial) 447 { 448 if (!t) t = new G4EnergyLossTablesHelper; 449 450 CPRWarning(); 451 if(aParticle != (const G4ParticleDefinition* 452 { 453 *t= GetTables(aParticle); 454 lastParticle = (G4ParticleDefinition*) aPa 455 oldIndex = -1 ; 456 } 457 const G4PhysicsTable* propertimeTable= t->th 458 if (!propertimeTable) { 459 ParticleHaveNoLoss(aParticle,"ProperTime") 460 return 0.0; 461 } 462 463 const G4double parlowen=0.4 , ppar=0.5-parlo 464 const G4double dToverT = 0.05 , facT = 1. -d 465 G4double timestart,timeend,deltatime,dTT; 466 G4bool isOut; 467 468 G4int materialIndex = (G4int)aMaterial->GetI 469 G4double scaledKineticEnergy = KineticEnergy 470 471 if (scaledKineticEnergy<t->theLowestKineticE 472 473 timestart = std::exp(ppar*std::log(scaled 474 (*propertimeTable)(materialInd 475 t->theLowestKineticEnergy,isOu 476 477 478 } else if (scaledKineticEnergy>t->theHighest 479 480 timestart = (*propertimeTable)(materialIn 481 t->theHighestKineticEnergy,isO 482 483 } else { 484 485 timestart = (*propertimeTable)(materialInd 486 scaledKineticEnergy,isOut); 487 488 } 489 490 dTT = (KineticEnergyStart - KineticEnergyEnd 491 492 if( dTT < dToverT ) 493 scaledKineticEnergy = facT*KineticEnergySt 494 else 495 scaledKineticEnergy = KineticEnergyEnd*t-> 496 497 if (scaledKineticEnergy<t->theLowestKineticE 498 499 timeend = std::exp(ppar*std::log(scaledKi 500 (*propertimeTable)(materialInd 501 t->theLowestKineticEnergy,isOu 502 503 504 } else if (scaledKineticEnergy>t->theHighest 505 506 timeend = (*propertimeTable)(materialInde 507 t->theHighestKineticEnergy,isO 508 509 } else { 510 511 timeend = (*propertimeTable)(materialIndex 512 scaledKineticEnergy,isOut); 513 514 } 515 516 deltatime = timestart - timeend ; 517 518 if( dTT < dToverT ) 519 deltatime *= dTT/dToverT ; 520 521 return deltatime/t->theMassRatio ; 522 } 523 524 //....oooOO0OOooo........oooOO0OOooo........oo 525 526 G4double G4EnergyLossTables::GetRange( 527 const G4ParticleDefinition *aParticle, 528 G4double KineticEnergy, 529 const G4Material *aMaterial) 530 { 531 if (!t) t = new G4EnergyLossTablesHelper; 532 533 CPRWarning(); 534 if(aParticle != (const G4ParticleDefinition* 535 { 536 *t= GetTables(aParticle); 537 lastParticle = (G4ParticleDefinition*) aPa 538 Chargesquare = (aParticle->GetPDGCharge()) 539 (aParticle->GetPDGCharge()) 540 QQPositron ; 541 oldIndex = -1 ; 542 } 543 const G4PhysicsTable* rangeTable= t->theRang 544 const G4PhysicsTable* dEdxTable= t->theDEDX 545 if (!rangeTable) { 546 ParticleHaveNoLoss(aParticle,"Range"); 547 return 0.0; 548 } 549 550 G4int materialIndex = (G4int)aMaterial->GetI 551 G4double scaledKineticEnergy = KineticEnergy 552 G4double Range; 553 G4bool isOut; 554 555 if (scaledKineticEnergy<t->theLowestKineticE 556 557 Range = std::sqrt(scaledKineticEnergy/t->t 558 (*rangeTable)(materialIndex)->GetV 559 t->theLowestKineticEnergy,isOut) 560 561 } else if (scaledKineticEnergy>t->theHighest 562 563 Range = (*rangeTable)(materialIndex)->GetV 564 t->theHighestKineticEnergy,isOut)+ 565 (scaledKineticEnergy-t->theHighest 566 (*dEdxTable)(materialIndex)->GetVa 567 t->theHighestKineticEnergy,isOut 568 569 } else { 570 571 Range = (*rangeTable)(materialIndex)->GetV 572 scaledKineticEnergy,isOut); 573 574 } 575 576 return Range/(Chargesquare*t->theMassRatio); 577 } 578 579 //....oooOO0OOooo........oooOO0OOooo........oo 580 581 G4double G4EnergyLossTables::GetPreciseEnergyF 582 const G4P 583 G4d 584 const G4M 585 // it returns the value of the kinetic energy 586 { 587 if (!t) t = new G4EnergyLossTablesHelper; 588 589 CPRWarning(); 590 if( aParticle != (const G4ParticleDefinition 591 { 592 *t= GetTables(aParticle); 593 lastParticle = (G4ParticleDefinition*) aPa 594 Chargesquare = (aParticle->GetPDGCharge()) 595 (aParticle->GetPDGCharge()) 596 QQPositron ; 597 oldIndex = -1 ; 598 } 599 const G4PhysicsTable* dEdxTable= t->theDEDX 600 const G4PhysicsTable* inverseRangeTable= t- 601 if (!inverseRangeTable) { 602 ParticleHaveNoLoss(aParticle,"InverseRange 603 return 0.0; 604 } 605 606 G4double scaledrange,scaledKineticEnergy ; 607 G4bool isOut ; 608 609 G4int materialIndex = (G4int)aMaterial->GetI 610 611 if(materialIndex != oldIndex) 612 { 613 oldIndex = materialIndex ; 614 rmin = (*inverseRangeTable)(materialIndex) 615 GetLowEdgeEnergy 616 rmax = (*inverseRangeTable)(materialIndex) 617 GetLowEdgeEnergy(t->theNumb 618 Thigh = (*inverseRangeTable)(materialIndex 619 GetValue(rmax,is 620 } 621 622 scaledrange = range*Chargesquare*t->theMassR 623 624 if(scaledrange < rmin) 625 { 626 scaledKineticEnergy = t->theLowestKineticE 627 scaledrange*scaledrange/(rm 628 } 629 else 630 { 631 if(scaledrange < rmax) 632 { 633 scaledKineticEnergy = (*inverseRangeTabl 634 GetValue( scaled 635 } 636 else 637 { 638 scaledKineticEnergy = Thigh + 639 (scaledrange-rmax)* 640 (*dEdxTable)(materialInd 641 GetValue(Thig 642 } 643 } 644 645 return scaledKineticEnergy/t->theMassRatio ; 646 } 647 648 //....oooOO0OOooo........oooOO0OOooo........oo 649 650 G4double G4EnergyLossTables::GetPreciseDEDX( 651 const G4ParticleDefinition *aParticle, 652 G4double KineticEnergy, 653 const G4Material *aMaterial) 654 { 655 if (!t) t = new G4EnergyLossTablesHelper; 656 657 CPRWarning(); 658 if( aParticle != (const G4ParticleDefinition 659 { 660 *t= GetTables(aParticle); 661 lastParticle = (G4ParticleDefinition*) aPa 662 Chargesquare = (aParticle->GetPDGCharge()) 663 (aParticle->GetPDGCharge()) 664 QQPositron ; 665 oldIndex = -1 ; 666 } 667 const G4PhysicsTable* dEdxTable= t->theDEDX 668 if (!dEdxTable) { 669 ParticleHaveNoLoss(aParticle,"dEdx"); 670 return 0.0; 671 } 672 673 G4int materialIndex = (G4int)aMaterial->GetI 674 G4double scaledKineticEnergy = KineticEnergy 675 G4double dEdx; 676 G4bool isOut; 677 678 if (scaledKineticEnergy<t->theLowestKineticE 679 680 dEdx = std::sqrt(scaledKineticEnergy/t->t 681 *(*dEdxTable)(materialIndex)->GetV 682 t->theLowestKineticEnergy,isOut) 683 684 } else if (scaledKineticEnergy>t->theHighest 685 686 dEdx = (*dEdxTable)(materialIndex)->GetVa 687 t->theHighestKineticEnergy,isOut); 688 689 } else { 690 691 dEdx = (*dEdxTable)(materialIndex)->GetV 692 scaledKineticEnergy, 693 694 } 695 696 return dEdx*Chargesquare; 697 } 698 699 //....oooOO0OOooo........oooOO0OOooo........oo 700 701 G4double G4EnergyLossTables::GetPreciseRangeF 702 const G4ParticleDefinition *aParticle, 703 G4double KineticEnergy, 704 const G4Material *aMaterial) 705 { 706 if (!t) t = new G4EnergyLossTablesHelper; 707 708 CPRWarning(); 709 if( aParticle != (const G4ParticleDefinition 710 { 711 *t= GetTables(aParticle); 712 lastParticle = (G4ParticleDefinition*) aPa 713 Chargesquare = (aParticle->GetPDGCharge()) 714 (aParticle->GetPDGCharge()) 715 QQPositron ; 716 oldIndex = -1 ; 717 } 718 const G4PhysicsTable* rangeTable= t->theRang 719 const G4PhysicsTable* dEdxTable= t->theDEDX 720 if (!rangeTable) { 721 ParticleHaveNoLoss(aParticle,"Range"); 722 return 0.0; 723 } 724 G4int materialIndex = (G4int)aMaterial->GetI 725 726 G4double Thighr = t->theHighestKineticEnergy 727 (*rangeTable)(materialIndex 728 GetLowEdgeEnergy(1) ; 729 730 G4double scaledKineticEnergy = KineticEnergy 731 G4double Range; 732 G4bool isOut; 733 734 if (scaledKineticEnergy<t->theLowestKineticE 735 736 Range = std::sqrt(scaledKineticEnergy/t->t 737 (*rangeTable)(materialIndex)->GetV 738 t->theLowestKineticEnergy,isOut) 739 740 } else if (scaledKineticEnergy>Thighr) { 741 742 Range = (*rangeTable)(materialIndex)->GetV 743 Thighr,isOut)+ 744 (scaledKineticEnergy-Thighr)/ 745 (*dEdxTable)(materialIndex)->GetVa 746 Thighr,isOut); 747 748 } else { 749 750 Range = (*rangeTable)(materialIndex)->Get 751 scaledKineticEnergy,isO 752 753 } 754 755 return Range/(Chargesquare*t->theMassRatio); 756 } 757 758 //....oooOO0OOooo........oooOO0OOooo........oo 759 760 G4double G4EnergyLossTables::GetDEDX( 761 const G4ParticleDefinition *aParticle, 762 G4double KineticEnergy, 763 const G4MaterialCutsCouple *couple, 764 G4bool check) 765 { 766 if (!t) t = new G4EnergyLossTablesHelper; 767 768 if(aParticle != (const G4ParticleDefinition* 769 { 770 *t= GetTables(aParticle); 771 lastParticle = (G4ParticleDefinition*) aPa 772 Chargesquare = (aParticle->GetPDGCharge()) 773 (aParticle->GetPDGCharge()) 774 QQPositron ; 775 oldIndex = -1 ; 776 } 777 const G4PhysicsTable* dEdxTable= t->theDEDX 778 779 if (!dEdxTable ) { 780 if (check) return G4LossTableManager::Inst 781 else ParticleHaveNoLoss(aParticle, " 782 return 0.0; 783 } 784 785 G4int materialIndex = couple->GetIndex(); 786 G4double scaledKineticEnergy = KineticEnergy 787 G4double dEdx; 788 G4bool isOut; 789 790 if (scaledKineticEnergy<t->theLowestKineticE 791 792 dEdx =(*dEdxTable)(materialIndex)->GetVal 793 t->theLowestKineticEnergy,isOut) 794 *std::sqrt(scaledKineticEnergy/t->t 795 796 } else if (scaledKineticEnergy>t->theHighest 797 798 dEdx = (*dEdxTable)(materialIndex)->GetVa 799 t->theHighestKineticEnergy,isOut); 800 801 } else { 802 803 dEdx = (*dEdxTable)(materialIndex)->GetVal 804 scaledKineticEnergy,isOut); 805 806 } 807 808 return dEdx*Chargesquare; 809 } 810 811 //....oooOO0OOooo........oooOO0OOooo........oo 812 813 G4double G4EnergyLossTables::GetRange( 814 const G4ParticleDefinition *aParticle, 815 G4double KineticEnergy, 816 const G4MaterialCutsCouple *couple, 817 G4bool check) 818 { 819 if (!t) t = new G4EnergyLossTablesHelper; 820 821 if(aParticle != (const G4ParticleDefinition* 822 { 823 *t= GetTables(aParticle); 824 lastParticle = (G4ParticleDefinition*) aPa 825 Chargesquare = (aParticle->GetPDGCharge()) 826 (aParticle->GetPDGCharge()) 827 QQPositron ; 828 oldIndex = -1 ; 829 } 830 const G4PhysicsTable* rangeTable= t->theRang 831 const G4PhysicsTable* dEdxTable= t->theDEDX 832 if (!rangeTable) { 833 if(check) return G4LossTableManager::Insta 834 else return DBL_MAX; 835 //ParticleHaveNoLoss(aParticle,"Range"); 836 } 837 838 G4int materialIndex = couple->GetIndex(); 839 G4double scaledKineticEnergy = KineticEnergy 840 G4double Range; 841 G4bool isOut; 842 843 if (scaledKineticEnergy<t->theLowestKineticE 844 845 Range = std::sqrt(scaledKineticEnergy/t->t 846 (*rangeTable)(materialIndex)->GetV 847 t->theLowestKineticEnergy,isOut) 848 849 } else if (scaledKineticEnergy>t->theHighest 850 851 Range = (*rangeTable)(materialIndex)->GetV 852 t->theHighestKineticEnergy,isOut)+ 853 (scaledKineticEnergy-t->theHighest 854 (*dEdxTable)(materialIndex)->GetVa 855 t->theHighestKineticEnergy,isOut 856 857 } else { 858 859 Range = (*rangeTable)(materialIndex)->GetV 860 scaledKineticEnergy,isOut); 861 862 } 863 864 return Range/(Chargesquare*t->theMassRatio); 865 } 866 867 //....oooOO0OOooo........oooOO0OOooo........oo 868 869 G4double G4EnergyLossTables::GetPreciseEnergyF 870 const G4P 871 G4d 872 const G4M 873 G4bool check) 874 // it returns the value of the kinetic energy 875 { 876 if (!t) t = new G4EnergyLossTablesHelper; 877 878 if( aParticle != (const G4ParticleDefinition 879 { 880 *t= GetTables(aParticle); 881 lastParticle = (G4ParticleDefinition*) aPa 882 Chargesquare = (aParticle->GetPDGCharge()) 883 (aParticle->GetPDGCharge()) 884 QQPositron ; 885 oldIndex = -1 ; 886 } 887 const G4PhysicsTable* dEdxTable= t->theDEDX 888 const G4PhysicsTable* inverseRangeTable= t- 889 890 if (!inverseRangeTable) { 891 if(check) return G4LossTableManager::Insta 892 else return DBL_MAX; 893 // else ParticleHaveNoLoss(aPartic 894 } 895 896 G4double scaledrange,scaledKineticEnergy ; 897 G4bool isOut ; 898 899 G4int materialIndex = couple->GetIndex() ; 900 901 if(materialIndex != oldIndex) 902 { 903 oldIndex = materialIndex ; 904 rmin = (*inverseRangeTable)(materialIndex) 905 GetLowEdgeEnergy 906 rmax = (*inverseRangeTable)(materialIndex) 907 GetLowEdgeEnergy(t->theNumb 908 Thigh = (*inverseRangeTable)(materialIndex 909 GetValue(rmax,is 910 } 911 912 scaledrange = range*Chargesquare*t->theMassR 913 914 if(scaledrange < rmin) 915 { 916 scaledKineticEnergy = t->theLowestKineticE 917 scaledrange*scaledrange/(rm 918 } 919 else 920 { 921 if(scaledrange < rmax) 922 { 923 scaledKineticEnergy = (*inverseRangeTabl 924 GetValue( scaled 925 } 926 else 927 { 928 scaledKineticEnergy = Thigh + 929 (scaledrange-rmax)* 930 (*dEdxTable)(materialInd 931 GetValue(Thig 932 } 933 } 934 935 return scaledKineticEnergy/t->theMassRatio ; 936 } 937 938 //....oooOO0OOooo........oooOO0OOooo........oo 939 940 G4double G4EnergyLossTables::GetPreciseDEDX( 941 const G4ParticleDefinition *aParticle, 942 G4double KineticEnergy, 943 const G4MaterialCutsCouple *couple) 944 { 945 if (!t) t = new G4EnergyLossTablesHelper; 946 947 if( aParticle != (const G4ParticleDefinition 948 { 949 *t= GetTables(aParticle); 950 lastParticle = (G4ParticleDefinition*) aPa 951 Chargesquare = (aParticle->GetPDGCharge()) 952 (aParticle->GetPDGCharge()) 953 QQPositron ; 954 oldIndex = -1 ; 955 } 956 const G4PhysicsTable* dEdxTable= t->theDEDX 957 if ( !dEdxTable ) 958 return G4LossTableManager::Instance()->Get 959 960 G4int materialIndex = couple->GetIndex(); 961 G4double scaledKineticEnergy = KineticEnergy 962 G4double dEdx; 963 G4bool isOut; 964 965 if (scaledKineticEnergy<t->theLowestKineticE 966 967 dEdx = std::sqrt(scaledKineticEnergy/t->t 968 *(*dEdxTable)(materialIndex)->GetV 969 t->theLowestKineticEnergy,isOut) 970 971 } else if (scaledKineticEnergy>t->theHighest 972 973 dEdx = (*dEdxTable)(materialIndex)->GetVa 974 t->theHighestKineticEnergy,isOut); 975 976 } else { 977 978 dEdx = (*dEdxTable)(materialIndex)->GetV 979 scaledKineticEnergy, 980 981 } 982 983 return dEdx*Chargesquare; 984 } 985 986 //....oooOO0OOooo........oooOO0OOooo........oo 987 988 G4double G4EnergyLossTables::GetPreciseRangeFr 989 const G4ParticleDefinition *aParticle, 990 G4double KineticEnergy, 991 const G4MaterialCutsCouple *couple) 992 { 993 if (!t) t = new G4EnergyLossTablesHelper; 994 995 if( aParticle != (const G4ParticleDefinition 996 { 997 *t= GetTables(aParticle); 998 lastParticle = (G4ParticleDefinition*) aPa 999 Chargesquare = (aParticle->GetPDGCharge()) 1000 (aParticle->GetPDGCharge() 1001 QQPositron ; 1002 oldIndex = -1 ; 1003 } 1004 const G4PhysicsTable* rangeTable= t->theRan 1005 const G4PhysicsTable* dEdxTable= t->theDED 1006 if ( !dEdxTable || !rangeTable) 1007 return G4LossTableManager::Instance()->Ge 1008 1009 G4int materialIndex = couple->GetIndex(); 1010 1011 G4double Thighr = t->theHighestKineticEnerg 1012 (*rangeTable)(materialInde 1013 GetLowEdgeEnergy(1) ; 1014 1015 G4double scaledKineticEnergy = KineticEnerg 1016 G4double Range; 1017 G4bool isOut; 1018 1019 if (scaledKineticEnergy<t->theLowestKinetic 1020 1021 Range = std::sqrt(scaledKineticEnergy/t-> 1022 (*rangeTable)(materialIndex)->Get 1023 t->theLowestKineticEnergy,isOut 1024 1025 } else if (scaledKineticEnergy>Thighr) { 1026 1027 Range = (*rangeTable)(materialIndex)->Get 1028 Thighr,isOut)+ 1029 (scaledKineticEnergy-Thighr)/ 1030 (*dEdxTable)(materialIndex)->GetV 1031 Thighr,isOut); 1032 1033 } else { 1034 1035 Range = (*rangeTable)(materialIndex)->Ge 1036 scaledKineticEnergy,is 1037 1038 } 1039 1040 return Range/(Chargesquare*t->theMassRatio) 1041 } 1042 1043 //....oooOO0OOooo........oooOO0OOooo........o 1044 1045 void G4EnergyLossTables::CPRWarning() 1046 { 1047 if (let_counter < let_max_num_warnings) { 1048 1049 G4cout << G4endl; 1050 G4cout << "##### G4EnergyLossTable WARNIN 1051 G4cout << "##### RESULTS ARE NOT GARANTEE 1052 G4cout << "##### Please, substitute G4Mat 1053 G4cout << "##### Obsolete interface will 1054 G4cout << G4endl; 1055 let_counter++; 1056 1057 } else if (let_counter == let_max_num_warni 1058 1059 G4cout << "##### G4EnergyLossTable WARNIN 1060 let_counter++; 1061 } 1062 } 1063 1064 //....oooOO0OOooo........oooOO0OOooo........o 1065 1066 void 1067 G4EnergyLossTables::ParticleHaveNoLoss(const 1068 const G4String& /*q*/) 1069 { 1070 } 1071 1072 //....oooOO0OOooo........oooOO0OOooo........o 1073