Geant4 Cross Reference |
>> 1 // This code implementation is the intellectual property of >> 2 // the GEANT4 collaboration. >> 3 // >> 4 // By copying, distributing or modifying the Program (or any work >> 5 // based on the Program) you indicate your acceptance of this statement, >> 6 // and all its terms. 1 // 7 // 2 // ******************************************* << 8 // $Id: G4EnergyLossTables.cc,v 1.12 2000/05/23 14:25:26 urban Exp $ 3 // * License and Disclaimer << 9 // GEANT4 tag $Name: geant4-02-00 $ 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 // 10 // 26 // << 11 // ------------------------------------------------------------------- 27 // << 28 // ------------------------------------------- << 29 // first version created by P.Urban , 06/04/19 12 // first version created by P.Urban , 06/04/1998 30 // modifications + "precise" functions added b 13 // modifications + "precise" functions added by L.Urban , 27/05/98 31 // modifications , TOF functions , 26/10/98, L 14 // modifications , TOF functions , 26/10/98, L.Urban 32 // cache mechanism in order to gain time, 11/0 15 // cache mechanism in order to gain time, 11/02/99, L.Urban 33 // bug fixed , 12/04/99 , L.Urban 16 // bug fixed , 12/04/99 , L.Urban 34 // 10.11.99: moved from RWT hash dictionary to << 17 // 10/11/99: moved from RWT hash dictionary to STL map, G.Barrand, M.Maire 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 // ------------------------------------------- 18 // ------------------------------------------------------------------- 42 19 43 #include "G4EnergyLossTables.hh" 20 #include "G4EnergyLossTables.hh" 44 #include "G4SystemOfUnits.hh" << 45 #include "G4MaterialCutsCouple.hh" << 46 #include "G4RegionStore.hh" << 47 #include "G4LossTableManager.hh" << 48 << 49 21 50 //....oooOO0OOooo........oooOO0OOooo........oo 22 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 51 23 52 G4EnergyLossTablesHelper *G4EnergyLossTables:: << 24 G4EnergyLossTablesHelper G4EnergyLossTables::t ; 53 G4EnergyLossTablesHelper *G4EnergyLossTables:: << 25 const G4ParticleDefinition* G4EnergyLossTables::lastParticle = NULL ; 54 G4ParticleDefinition* G4EnergyLossTables::last << 26 G4double G4EnergyLossTables::QQPositron = eplus*eplus ; 55 G4double G4EnergyLossTables::QQPositron = 1.0; << 56 G4double G4EnergyLossTables::Chargesquare ; 27 G4double G4EnergyLossTables::Chargesquare ; 57 G4int G4EnergyLossTables::oldIndex = -1 ; 28 G4int G4EnergyLossTables::oldIndex = -1 ; 58 G4double G4EnergyLossTables::rmin = 0. ; 29 G4double G4EnergyLossTables::rmin = 0. ; 59 G4double G4EnergyLossTables::rmax = 0. ; 30 G4double G4EnergyLossTables::rmax = 0. ; 60 G4double G4EnergyLossTables::Thigh = 0. ; 31 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 32 69 G4EnergyLossTablesHelper::G4EnergyLossTablesHe << 33 G4EnergyLossTables::helper_map G4EnergyLossTables::dict; 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 34 102 //....oooOO0OOooo........oooOO0OOooo........oo 35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 103 36 104 void G4EnergyLossTables::Register( 37 void G4EnergyLossTables::Register( 105 const G4ParticleDefinition* p, 38 const G4ParticleDefinition* p, 106 const G4PhysicsTable* tDEDX, 39 const G4PhysicsTable* tDEDX, 107 const G4PhysicsTable* tRange, 40 const G4PhysicsTable* tRange, 108 const G4PhysicsTable* tInverseRange, 41 const G4PhysicsTable* tInverseRange, 109 const G4PhysicsTable* tLabTime, 42 const G4PhysicsTable* tLabTime, 110 const G4PhysicsTable* tProperTime, 43 const G4PhysicsTable* tProperTime, 111 G4double lowestKineticEnergy, 44 G4double lowestKineticEnergy, 112 G4double highestKineticEnergy, 45 G4double highestKineticEnergy, 113 G4double massRatio, 46 G4double massRatio, 114 G4int NumberOfBins) 47 G4int NumberOfBins) 115 { 48 { 116 if (!dict) dict = new G4EnergyLossTables::he << 49 dict[p]= G4EnergyLossTablesHelper(tDEDX, tRange,tInverseRange, 117 if (!null_loss) null_loss = new G4EnergyLoss << 118 if (!t) t = new G4EnergyLossTablesHelper; << 119 << 120 (*dict)[p]= G4EnergyLossTablesHelper(tDEDX, << 121 tLabTime,tProperTime,lowes 50 tLabTime,tProperTime,lowestKineticEnergy, 122 highestKineticEnergy, massRatio,Number 51 highestKineticEnergy, massRatio,NumberOfBins); 123 << 52 124 *t = GetTables(p) ; // important for cach << 53 t = GetTables(p) ; // important for cache !!!!! 125 lastParticle = (G4ParticleDefinition*) p ; << 54 lastParticle = p ; 126 Chargesquare = (p->GetPDGCharge())*(p->GetPD 55 Chargesquare = (p->GetPDGCharge())*(p->GetPDGCharge())/ 127 QQPositron ; 56 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 } 57 } 647 58 648 //....oooOO0OOooo........oooOO0OOooo........oo 59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 649 60 650 G4double G4EnergyLossTables::GetPreciseDEDX( 61 G4double G4EnergyLossTables::GetPreciseDEDX( 651 const G4ParticleDefinition *aParticle, 62 const G4ParticleDefinition *aParticle, 652 G4double KineticEnergy, 63 G4double KineticEnergy, 653 const G4Material *aMaterial) << 64 G4Material *aMaterial) 654 { 65 { 655 if (!t) t = new G4EnergyLossTablesHelper; << 66 if( aParticle != lastParticle) 656 << 657 CPRWarning(); << 658 if( aParticle != (const G4ParticleDefinition << 659 { 67 { 660 *t= GetTables(aParticle); << 68 t= GetTables(aParticle); 661 lastParticle = (G4ParticleDefinition*) aPa << 69 lastParticle = aParticle; 662 Chargesquare = (aParticle->GetPDGCharge()) 70 Chargesquare = (aParticle->GetPDGCharge())* 663 (aParticle->GetPDGCharge()) 71 (aParticle->GetPDGCharge())/ 664 QQPositron ; 72 QQPositron ; 665 oldIndex = -1 ; << 666 } << 667 const G4PhysicsTable* dEdxTable= t->theDEDX << 668 if (!dEdxTable) { << 669 ParticleHaveNoLoss(aParticle,"dEdx"); << 670 return 0.0; << 671 } 73 } >> 74 const G4PhysicsTable* dEdxTable= t.theDEDXTable; 672 75 673 G4int materialIndex = (G4int)aMaterial->GetI << 76 G4int materialIndex = aMaterial->GetIndex(); 674 G4double scaledKineticEnergy = KineticEnergy << 77 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio; 675 G4double dEdx; 78 G4double dEdx; 676 G4bool isOut; 79 G4bool isOut; 677 80 678 if (scaledKineticEnergy<t->theLowestKineticE << 81 if (scaledKineticEnergy<t.theLowestKineticEnergy) { 679 82 680 dEdx = std::sqrt(scaledKineticEnergy/t->t << 83 dEdx = sqrt(scaledKineticEnergy/t.theLowestKineticEnergy) 681 *(*dEdxTable)(materialIndex)->GetV 84 *(*dEdxTable)(materialIndex)->GetValue( 682 t->theLowestKineticEnergy,isOut) << 85 t.theLowestKineticEnergy,isOut); 683 86 684 } else if (scaledKineticEnergy>t->theHighest << 87 } else if (scaledKineticEnergy>t.theHighestKineticEnergy) { 685 88 686 dEdx = (*dEdxTable)(materialIndex)->GetVa 89 dEdx = (*dEdxTable)(materialIndex)->GetValue( 687 t->theHighestKineticEnergy,isOut); << 90 t.theHighestKineticEnergy,isOut); 688 91 689 } else { 92 } else { 690 << 93 691 dEdx = (*dEdxTable)(materialIndex)->GetV 94 dEdx = (*dEdxTable)(materialIndex)->GetValue( 692 scaledKineticEnergy, 95 scaledKineticEnergy,isOut) ; 693 96 694 } 97 } 695 98 696 return dEdx*Chargesquare; 99 return dEdx*Chargesquare; 697 } 100 } 698 101 699 //....oooOO0OOooo........oooOO0OOooo........oo 102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 700 103 701 G4double G4EnergyLossTables::GetPreciseRangeF 104 G4double G4EnergyLossTables::GetPreciseRangeFromEnergy( 702 const G4ParticleDefinition *aParticle, 105 const G4ParticleDefinition *aParticle, 703 G4double KineticEnergy, 106 G4double KineticEnergy, 704 const G4Material *aMaterial) << 107 G4Material *aMaterial) 705 { 108 { 706 if (!t) t = new G4EnergyLossTablesHelper; << 109 if( aParticle != lastParticle) 707 << 708 CPRWarning(); << 709 if( aParticle != (const G4ParticleDefinition << 710 { 110 { 711 *t= GetTables(aParticle); << 111 t= GetTables(aParticle); 712 lastParticle = (G4ParticleDefinition*) aPa << 112 lastParticle = aParticle; 713 Chargesquare = (aParticle->GetPDGCharge()) 113 Chargesquare = (aParticle->GetPDGCharge())* 714 (aParticle->GetPDGCharge()) 114 (aParticle->GetPDGCharge())/ 715 QQPositron ; 115 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 } 116 } >> 117 const G4PhysicsTable* rangeTable= t.theRangeTable; >> 118 const G4PhysicsTable* dEdxTable= t.theDEDXTable; 754 119 755 return Range/(Chargesquare*t->theMassRatio); << 120 G4int materialIndex = aMaterial->GetIndex(); 756 } << 757 121 758 //....oooOO0OOooo........oooOO0OOooo........oo << 122 G4double Thighr = t.theHighestKineticEnergy*t.theLowestKineticEnergy/ 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 123 (*rangeTable)(materialIndex)-> 1013 GetLowEdgeEnergy(1) ; 124 GetLowEdgeEnergy(1) ; 1014 125 1015 G4double scaledKineticEnergy = KineticEnerg << 126 G4double scaledKineticEnergy = KineticEnergy*t.theMassRatio; 1016 G4double Range; 127 G4double Range; 1017 G4bool isOut; 128 G4bool isOut; 1018 129 1019 if (scaledKineticEnergy<t->theLowestKinetic << 130 if (scaledKineticEnergy<t.theLowestKineticEnergy) { 1020 131 1021 Range = std::sqrt(scaledKineticEnergy/t-> << 132 Range = sqrt(scaledKineticEnergy/t.theLowestKineticEnergy)* 1022 (*rangeTable)(materialIndex)->Get 133 (*rangeTable)(materialIndex)->GetValue( 1023 t->theLowestKineticEnergy,isOut << 134 t.theLowestKineticEnergy,isOut); 1024 135 1025 } else if (scaledKineticEnergy>Thighr) { 136 } else if (scaledKineticEnergy>Thighr) { 1026 137 1027 Range = (*rangeTable)(materialIndex)->Get 138 Range = (*rangeTable)(materialIndex)->GetValue( 1028 Thighr,isOut)+ 139 Thighr,isOut)+ 1029 (scaledKineticEnergy-Thighr)/ 140 (scaledKineticEnergy-Thighr)/ 1030 (*dEdxTable)(materialIndex)->GetV 141 (*dEdxTable)(materialIndex)->GetValue( 1031 Thighr,isOut); 142 Thighr,isOut); 1032 143 1033 } else { 144 } else { 1034 << 145 1035 Range = (*rangeTable)(materialIndex)->Ge 146 Range = (*rangeTable)(materialIndex)->GetValue( 1036 scaledKineticEnergy,is 147 scaledKineticEnergy,isOut) ; 1037 148 1038 } 149 } 1039 150 1040 return Range/(Chargesquare*t->theMassRatio) << 151 return Range/(Chargesquare*t.theMassRatio); 1041 } 152 } 1042 153 1043 //....oooOO0OOooo........oooOO0OOooo........o 154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 1044 155 1045 void G4EnergyLossTables::CPRWarning() << 1046 { << 1047 if (let_counter < let_max_num_warnings) { << 1048 156 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 157 1057 } else if (let_counter == let_max_num_warni << 1058 << 1059 G4cout << "##### G4EnergyLossTable WARNIN << 1060 let_counter++; << 1061 } << 1062 } << 1063 158 1064 //....oooOO0OOooo........oooOO0OOooo........o << 1065 << 1066 void << 1067 G4EnergyLossTables::ParticleHaveNoLoss(const << 1068 const G4String& /*q*/) << 1069 { << 1070 } << 1071 << 1072 //....oooOO0OOooo........oooOO0OOooo........o << 1073 159