Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * * 21 // * Parts of this code which have been developed by QinetiQ Ltd * 22 // * under contract to the European Space Agency (ESA) are the * 23 // * intellectual property of ESA. Rights to use, copy, modify and * 24 // * redistribute this software for general public use are granted * 25 // * in compliance with any licensing, distribution and development * 26 // * policy adopted by the Geant4 Collaboration. This code has been * 27 // * written by QinetiQ Ltd for the European Space Agency, under ESA * 28 // * contract 17191/03/NL/LvH (Aurora Programme). * 29 // * * 30 // * By using, copying, modifying or distributing the software (or * 31 // * any work based on the software) you agree to acknowledge its * 32 // * use in resulting scientific publications, and indicate your * 33 // * acceptance of all terms of the Geant4 Software license. * 34 // ******************************************************************** 35 // 36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 37 // 38 // MODULE: G4WilsonAblationModel.cc 39 // 40 // Version: 1.0 41 // Date: 08/12/2009 42 // Author: P R Truscott 43 // Organisation: QinetiQ Ltd, UK 44 // Customer: ESA/ESTEC, NOORDWIJK 45 // Contract: 17191/03/NL/LvH 46 // 47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 48 // 49 // CHANGE HISTORY 50 // -------------- 51 // 52 // 6 October 2003, P R Truscott, QinetiQ Ltd, UK 53 // Created. 54 // 55 // 15 March 2004, P R Truscott, QinetiQ Ltd, UK 56 // Beta release 57 // 58 // 08 December 2009, P R Truscott, QinetiQ Ltd, UK 59 // Ver 1.0 60 // Updated as a result of changes in the G4Evaporation classes. These changes 61 // affect mostly SelectSecondariesByEvaporation, and now you have variables 62 // associated with the evaporation model which can be changed: 63 // OPTxs to select the inverse cross-section 64 // OPTxs = 0 => Dostrovski's parameterization 65 // OPTxs = 1 or 2 => Chatterjee's paramaterization 66 // OPTxs = 3 or 4 => Kalbach's parameterization 67 // useSICB => use superimposed Coulomb Barrier for inverse cross 68 // sections 69 // Other problem found with G4Fragment definition using Lorentz vector and 70 // **G4ParticleDefinition**. This does not allow A and Z to be defined for the 71 // fragment for some reason. Now the fragment is defined more explicitly: 72 // G4Fragment *fragment = new G4Fragment(A, Z, lorentzVector); 73 // to avoid this quirk. Bug found in SelectSecondariesByDefault: *type is now 74 // equated to evapType[i] whereas previously it was equated to fragType[i]. 75 // 76 // 06 August 2015, A. Ribon, CERN 77 // Migrated std::exp and std::pow to the faster G4Exp and G4Pow. 78 // 79 // 09 June 2017, C. Mancini Terracciano, INFN 80 // Fixed bug on the initialization of Photon Evaporation model 81 // 82 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 83 //////////////////////////////////////////////////////////////////////////////// 84 // 85 #include <iomanip> 86 #include <numeric> 87 88 #include "G4WilsonAblationModel.hh" 89 #include "G4PhysicalConstants.hh" 90 #include "G4SystemOfUnits.hh" 91 #include "Randomize.hh" 92 #include "G4ParticleTable.hh" 93 #include "G4IonTable.hh" 94 #include "G4Alpha.hh" 95 #include "G4He3.hh" 96 #include "G4Triton.hh" 97 #include "G4Deuteron.hh" 98 #include "G4Proton.hh" 99 #include "G4Neutron.hh" 100 #include "G4AlphaEvaporationChannel.hh" 101 #include "G4He3EvaporationChannel.hh" 102 #include "G4TritonEvaporationChannel.hh" 103 #include "G4DeuteronEvaporationChannel.hh" 104 #include "G4ProtonEvaporationChannel.hh" 105 #include "G4NeutronEvaporationChannel.hh" 106 #include "G4PhotonEvaporation.hh" 107 #include "G4LorentzVector.hh" 108 #include "G4VEvaporationChannel.hh" 109 110 #include "G4Exp.hh" 111 #include "G4Pow.hh" 112 113 #include "G4PhysicsModelCatalog.hh" 114 115 //////////////////////////////////////////////////////////////////////////////// 116 // 117 G4WilsonAblationModel::G4WilsonAblationModel() 118 { 119 // 120 // 121 // Send message to stdout to advise that the G4Abrasion model is being used. 122 // 123 PrintWelcomeMessage(); 124 // 125 // 126 // Set the default verbose level to 0 - no output. 127 // 128 verboseLevel = 0; 129 // 130 // 131 // Set the binding energy per nucleon .... did I mention that this is a crude 132 // model for nuclear de-excitation? 133 // 134 B = 10.0 * MeV; 135 // 136 // 137 // It is possuble to switch off secondary particle production (other than the 138 // final nuclear fragment). The default is on. 139 // 140 produceSecondaries = true; 141 // 142 // 143 // Now we need to define the decay modes. We're using the G4Evaporation model 144 // to help determine the kinematics of the decay. 145 // 146 nFragTypes = 6; 147 fragType[0] = G4Alpha::Alpha(); 148 fragType[1] = G4He3::He3(); 149 fragType[2] = G4Triton::Triton(); 150 fragType[3] = G4Deuteron::Deuteron(); 151 fragType[4] = G4Proton::Proton(); 152 fragType[5] = G4Neutron::Neutron(); 153 for(G4int i=0; i<200; ++i) { fSig[i] = 0.0; } 154 // 155 // 156 // Set verboseLevel default to no output. 157 // 158 verboseLevel = 0; 159 theChannelFactory = new G4EvaporationFactory(new G4PhotonEvaporation()); 160 theChannels = theChannelFactory->GetChannel(); 161 // 162 // 163 // Set defaults for evaporation classes. These can be overridden by user 164 // "set" methods. 165 // 166 OPTxs = 3; 167 useSICB = false; 168 fragmentVector = 0; 169 170 secID = G4PhysicsModelCatalog::GetModelID("model_G4WilsonAblationModel"); 171 } 172 //////////////////////////////////////////////////////////////////////////////// 173 // 174 G4WilsonAblationModel::~G4WilsonAblationModel() 175 {} 176 177 //////////////////////////////////////////////////////////////////////////////// 178 // 179 G4FragmentVector *G4WilsonAblationModel::BreakItUp 180 (const G4Fragment &theNucleus) 181 { 182 // 183 // 184 // Initilise the pointer to the G4FragmentVector used to return the information 185 // about the breakup. 186 // 187 fragmentVector = new G4FragmentVector; 188 fragmentVector->clear(); 189 // 190 // 191 // Get the A, Z and excitation of the nucleus. 192 // 193 G4int A = theNucleus.GetA_asInt(); 194 G4int Z = theNucleus.GetZ_asInt(); 195 G4double ex = theNucleus.GetExcitationEnergy(); 196 if (verboseLevel >= 2) 197 { 198 G4cout <<"oooooooooooooooooooooooooooooooooooooooo" 199 <<"oooooooooooooooooooooooooooooooooooooooo" 200 <<G4endl; 201 G4cout.precision(6); 202 G4cout <<"IN G4WilsonAblationModel" <<G4endl; 203 G4cout <<"Initial prefragment A=" <<A 204 <<", Z=" <<Z 205 <<", excitation energy = " <<ex/MeV <<" MeV" 206 <<G4endl; 207 } 208 // 209 // 210 // Check that there is a nucleus to speak of. It's possible there isn't one 211 // or its just a proton or neutron. In either case, the excitation energy 212 // (from the Lorentz vector) is not used. 213 // 214 if (A == 0) 215 { 216 if (verboseLevel >= 2) 217 { 218 G4cout <<"No nucleus to decay" <<G4endl; 219 G4cout <<"oooooooooooooooooooooooooooooooooooooooo" 220 <<"oooooooooooooooooooooooooooooooooooooooo" 221 <<G4endl; 222 } 223 return fragmentVector; 224 } 225 else if (A == 1) 226 { 227 G4LorentzVector lorentzVector = theNucleus.GetMomentum(); 228 lorentzVector.setE(lorentzVector.e()-ex+10.0*eV); 229 if (Z == 0) 230 { 231 G4Fragment *fragment = new G4Fragment(lorentzVector,G4Neutron::Neutron()); 232 if (fragment != nullptr) { fragment->SetCreatorModelID(secID); } 233 fragmentVector->push_back(fragment); 234 } 235 else 236 { 237 G4Fragment *fragment = new G4Fragment(lorentzVector,G4Proton::Proton()); 238 if (fragment != nullptr) { fragment->SetCreatorModelID(secID); } 239 fragmentVector->push_back(fragment); 240 } 241 if (verboseLevel >= 2) 242 { 243 G4cout <<"Final fragment is in fact only a nucleon) :" <<G4endl; 244 G4cout <<(*fragmentVector)[0] <<G4endl; 245 G4cout <<"oooooooooooooooooooooooooooooooooooooooo" 246 <<"oooooooooooooooooooooooooooooooooooooooo" 247 <<G4endl; 248 } 249 return fragmentVector; 250 } 251 // 252 // 253 // Then the number of nucleons ablated (either as nucleons or light nuclear 254 // fragments) is based on a simple argument for the binding energy per nucleon. 255 // 256 G4int DAabl = (G4int) (ex / B); 257 if (DAabl > A) DAabl = A; 258 // The following lines are no longer accurate given we now treat the final fragment 259 // if (verboseLevel >= 2) 260 // G4cout <<"Number of nucleons ejected = " <<DAabl <<G4endl; 261 262 // 263 // 264 // Determine the nuclear fragment from the ablation process by sampling the 265 // Rudstam equation. 266 // 267 G4int AF = A - DAabl; 268 G4int ZF = 0; 269 270 if (AF > 0) 271 { 272 G4Pow* g4calc = G4Pow::GetInstance(); 273 G4double AFd = (G4double) AF; 274 G4double R = 11.8 / g4calc->powZ(AF, 0.45); 275 G4int minZ = std::max(1, Z - DAabl); 276 // 277 // 278 // Here we define an integral probability distribution based on the Rudstam 279 // equation assuming a constant AF. 280 // 281 G4int zmax = std::min(199, Z); 282 G4double sum = 0.0; 283 for (ZF=minZ; ZF<=zmax; ++ZF) 284 { 285 sum += G4Exp(-R*g4calc->powA(std::abs(ZF - 0.486*AFd + 3.8E-04*AFd*AFd),1.5)); 286 fSig[ZF] = sum; 287 } 288 // 289 // 290 // Now sample that distribution to determine a value for ZF. 291 // 292 sum *= G4UniformRand(); 293 for (ZF=minZ; ZF<=zmax; ++ZF) { 294 if(sum <= fSig[ZF]) { break; } 295 } 296 } 297 G4int DZabl = Z - ZF; 298 // 299 // 300 // Now determine the nucleons or nuclei which have bee ablated. The preference 301 // is for the production of alphas, then other nuclei in order of decreasing 302 // binding energy. The energies assigned to the products of the decay are 303 // provisional for the moment (the 10eV is just to avoid errors with negative 304 // excitation energies due to rounding). 305 // 306 G4double totalEpost = 0.0; 307 evapType.clear(); 308 for (G4int ift=0; ift<nFragTypes; ift++) 309 { 310 G4ParticleDefinition *type = fragType[ift]; 311 G4double n = std::floor((G4double) DAabl / type->GetBaryonNumber() + 1.0E-10); 312 G4double n1 = 1.0E+10; 313 if (fragType[ift]->GetPDGCharge() > 0.0) 314 n1 = std::floor((G4double) DZabl / type->GetPDGCharge() + 1.0E-10); 315 if (n > n1) n = n1; 316 if (n > 0.0) 317 { 318 G4double mass = type->GetPDGMass(); 319 for (G4int j=0; j<(G4int) n; j++) 320 { 321 totalEpost += mass; 322 evapType.push_back(type); 323 } 324 DAabl -= (G4int) (n * type->GetBaryonNumber() + 1.0E-10); 325 DZabl -= (G4int) (n * type->GetPDGCharge() + 1.0E-10); 326 } 327 } 328 // 329 // 330 // Determine the properties of the final nuclear fragment. Note that if 331 // the final fragment is predicted to have a nucleon number of zero, then 332 // really it's the particle last in the vector evapType which becomes the 333 // final fragment. Therefore delete this from the vector if this is the 334 // case. 335 // 336 G4double massFinalFrag = 0.0; 337 if (AF > 0) 338 massFinalFrag = G4ParticleTable::GetParticleTable()->GetIonTable()-> 339 GetIonMass(ZF,AF); 340 else 341 { 342 G4ParticleDefinition *type = evapType[evapType.size()-1]; 343 AF = type->GetBaryonNumber(); 344 ZF = (G4int) (type->GetPDGCharge() + 1.0E-10); 345 evapType.erase(evapType.end()-1); 346 } 347 totalEpost += massFinalFrag; 348 // 349 // 350 // Provide verbose output on the nuclear fragment if requested. 351 // 352 if (verboseLevel >= 2) 353 { 354 G4cout <<"Final fragment A=" <<AF 355 <<", Z=" <<ZF 356 <<G4endl; 357 for (G4int ift=0; ift<nFragTypes; ift++) 358 { 359 G4ParticleDefinition *type = fragType[ift]; 360 G4long n = std::count(evapType.cbegin(),evapType.cend(),type); 361 if (n > 0) 362 G4cout <<"Particle type: " <<std::setw(10) <<type->GetParticleName() 363 <<", number of particles emitted = " <<n <<G4endl; 364 } 365 } 366 // 367 // Add the total energy from the fragment. Note that the fragment is assumed 368 // to be de-excited and does not undergo photo-evaporation .... I did mention 369 // this is a bit of a crude model? 370 // 371 G4double massPreFrag = theNucleus.GetGroundStateMass(); 372 G4double totalEpre = massPreFrag + ex; 373 G4double excess = totalEpre - totalEpost; 374 // G4Fragment *resultNucleus(theNucleus); 375 G4Fragment *resultNucleus = new G4Fragment(A, Z, theNucleus.GetMomentum()); 376 G4ThreeVector boost(0.0,0.0,0.0); 377 std::size_t nEvap = 0; 378 if (produceSecondaries && evapType.size()>0) 379 { 380 if (excess > 0.0) 381 { 382 SelectSecondariesByEvaporation (resultNucleus); 383 nEvap = fragmentVector->size(); 384 boost = resultNucleus->GetMomentum().findBoostToCM(); 385 if (evapType.size() > 0) 386 SelectSecondariesByDefault (boost); 387 } 388 else 389 SelectSecondariesByDefault(G4ThreeVector(0.0,0.0,0.0)); 390 } 391 392 if (AF > 0) 393 { 394 G4double mass = G4ParticleTable::GetParticleTable()->GetIonTable()-> 395 GetIonMass(ZF,AF); 396 G4double e = mass + 10.0*eV; 397 G4double p = std::sqrt(e*e-mass*mass); 398 G4ThreeVector direction(0.0,0.0,1.0); 399 G4LorentzVector lorentzVector = G4LorentzVector(direction*p, e); 400 lorentzVector.boost(-boost); 401 G4Fragment* frag = new G4Fragment(AF, ZF, lorentzVector); 402 if (frag != nullptr) { frag->SetCreatorModelID(secID); } 403 fragmentVector->push_back(frag); 404 } 405 delete resultNucleus; 406 // 407 // 408 // Provide verbose output on the ablation products if requested. 409 // 410 if (verboseLevel >= 2) 411 { 412 if (nEvap > 0) 413 { 414 G4cout <<"----------------------" <<G4endl; 415 G4cout <<"Evaporated particles :" <<G4endl; 416 G4cout <<"----------------------" <<G4endl; 417 } 418 std::size_t ie = 0; 419 for (auto iter = fragmentVector->cbegin(); 420 iter != fragmentVector->cend(); ++iter) 421 { 422 if (ie == nEvap) 423 { 424 // G4cout <<*iter <<G4endl; 425 G4cout <<"---------------------------------" <<G4endl; 426 G4cout <<"Particles from default emission :" <<G4endl; 427 G4cout <<"---------------------------------" <<G4endl; 428 } 429 G4cout <<*iter <<G4endl; 430 } 431 G4cout <<"oooooooooooooooooooooooooooooooooooooooo" 432 <<"oooooooooooooooooooooooooooooooooooooooo" 433 <<G4endl; 434 } 435 436 return fragmentVector; 437 } 438 //////////////////////////////////////////////////////////////////////////////// 439 // 440 void G4WilsonAblationModel::SelectSecondariesByEvaporation 441 (G4Fragment *intermediateNucleus) 442 { 443 G4Fragment theResidualNucleus = *intermediateNucleus; 444 G4bool evaporate = true; 445 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko 446 while (evaporate && evapType.size() != 0) 447 { 448 // 449 // 450 // Here's the cheaky bit. We're hijacking the G4Evaporation model, in order to 451 // more accurately sample to kinematics, but the species of the nuclear 452 // fragments will be the ones of our choosing as above. 453 // 454 std::vector <G4VEvaporationChannel*> theChannels1; 455 theChannels1.clear(); 456 std::vector <G4VEvaporationChannel*>::iterator i; 457 VectorOfFragmentTypes::iterator iter; 458 std::vector <VectorOfFragmentTypes::iterator> iters; 459 iters.clear(); 460 iter = std::find(evapType.begin(), evapType.end(), G4Alpha::Alpha()); 461 if (iter != evapType.end()) 462 { 463 theChannels1.push_back(new G4AlphaEvaporationChannel); 464 i = theChannels1.end() - 1; 465 (*i)->SetOPTxs(OPTxs); 466 (*i)->UseSICB(useSICB); 467 // (*i)->Initialize(theResidualNucleus); 468 iters.push_back(iter); 469 } 470 iter = std::find(evapType.begin(), evapType.end(), G4He3::He3()); 471 if (iter != evapType.end()) 472 { 473 theChannels1.push_back(new G4He3EvaporationChannel); 474 i = theChannels1.end() - 1; 475 (*i)->SetOPTxs(OPTxs); 476 (*i)->UseSICB(useSICB); 477 // (*i)->Initialize(theResidualNucleus); 478 iters.push_back(iter); 479 } 480 iter = std::find(evapType.begin(), evapType.end(), G4Triton::Triton()); 481 if (iter != evapType.end()) 482 { 483 theChannels1.push_back(new G4TritonEvaporationChannel); 484 i = theChannels1.end() - 1; 485 (*i)->SetOPTxs(OPTxs); 486 (*i)->UseSICB(useSICB); 487 // (*i)->Initialize(theResidualNucleus); 488 iters.push_back(iter); 489 } 490 iter = std::find(evapType.begin(), evapType.end(), G4Deuteron::Deuteron()); 491 if (iter != evapType.end()) 492 { 493 theChannels1.push_back(new G4DeuteronEvaporationChannel); 494 i = theChannels1.end() - 1; 495 (*i)->SetOPTxs(OPTxs); 496 (*i)->UseSICB(useSICB); 497 // (*i)->Initialize(theResidualNucleus); 498 iters.push_back(iter); 499 } 500 iter = std::find(evapType.begin(), evapType.end(), G4Proton::Proton()); 501 if (iter != evapType.end()) 502 { 503 theChannels1.push_back(new G4ProtonEvaporationChannel); 504 i = theChannels1.end() - 1; 505 (*i)->SetOPTxs(OPTxs); 506 (*i)->UseSICB(useSICB); 507 // (*i)->Initialize(theResidualNucleus); 508 iters.push_back(iter); 509 } 510 iter = std::find(evapType.begin(), evapType.end(), G4Neutron::Neutron()); 511 if (iter != evapType.end()) 512 { 513 theChannels1.push_back(new G4NeutronEvaporationChannel); 514 i = theChannels1.end() - 1; 515 (*i)->SetOPTxs(OPTxs); 516 (*i)->UseSICB(useSICB); 517 // (*i)->Initialize(theResidualNucleus); 518 iters.push_back(iter); 519 } 520 std::size_t nChannels = theChannels1.size(); 521 522 G4double totalProb = 0.0; 523 G4int ich = 0; 524 G4double probEvapType[6] = {0.0}; 525 for (auto iterEv=theChannels1.cbegin(); 526 iterEv!=theChannels1.cend(); ++iterEv) { 527 totalProb += (*iterEv)->GetEmissionProbability(intermediateNucleus); 528 probEvapType[ich] = totalProb; 529 ++ich; 530 } 531 if (totalProb > 0.0) { 532 // 533 // 534 // The emission probability for at least one of the evaporation channels is 535 // positive, therefore work out which one should be selected and decay 536 // the nucleus. 537 // 538 G4double xi = totalProb*G4UniformRand(); 539 std::size_t ii = 0; 540 for (ii=0; ii<nChannels; ++ii) 541 { 542 if (xi < probEvapType[ii]) { break; } 543 } 544 if (ii >= nChannels) { ii = nChannels - 1; } 545 G4FragmentVector *evaporationResult = theChannels1[ii]-> 546 BreakUpFragment(intermediateNucleus); 547 if ((*evaporationResult)[0] != nullptr) 548 { 549 (*evaporationResult)[0]->SetCreatorModelID(secID); 550 } 551 fragmentVector->push_back((*evaporationResult)[0]); 552 intermediateNucleus = (*evaporationResult)[1]; 553 delete evaporationResult; 554 } 555 else 556 { 557 // 558 // 559 // Probability for further evaporation is nil so have to escape from this 560 // routine and set the energies of the secondaries to 10eV. 561 // 562 evaporate = false; 563 } 564 } 565 566 return; 567 } 568 //////////////////////////////////////////////////////////////////////////////// 569 // 570 void G4WilsonAblationModel::SelectSecondariesByDefault (G4ThreeVector boost) 571 { 572 for (std::size_t i=0; i<evapType.size(); ++i) 573 { 574 G4ParticleDefinition *type = evapType[i]; 575 G4double mass = type->GetPDGMass(); 576 G4double e = mass + 10.0*eV; 577 G4double p = std::sqrt(e*e-mass*mass); 578 G4double costheta = 2.0*G4UniformRand() - 1.0; 579 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta)); 580 G4double phi = twopi * G4UniformRand() * rad; 581 G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta); 582 G4LorentzVector lorentzVector = G4LorentzVector(direction*p, e); 583 lorentzVector.boost(-boost); 584 // Possibility that the following line is not correctly carrying over A and Z 585 // from particle definition. Force values. PRT 03/12/2009. 586 // G4Fragment *fragment = 587 // new G4Fragment(lorentzVector, type); 588 G4int A = type->GetBaryonNumber(); 589 G4int Z = (G4int) (type->GetPDGCharge() + 1.0E-10); 590 G4Fragment *fragment = 591 new G4Fragment(A, Z, lorentzVector); 592 if (fragment != nullptr) { fragment->SetCreatorModelID(secID); } 593 fragmentVector->push_back(fragment); 594 } 595 } 596 //////////////////////////////////////////////////////////////////////////////// 597 // 598 void G4WilsonAblationModel::PrintWelcomeMessage () 599 { 600 G4cout <<G4endl; 601 G4cout <<" *****************************************************************" 602 <<G4endl; 603 G4cout <<" Nuclear ablation model for nuclear-nuclear interactions activated" 604 <<G4endl; 605 G4cout <<" (Written by QinetiQ Ltd for the European Space Agency)" 606 <<G4endl; 607 G4cout <<" !!! WARNING: This model is not well validation and should not be used for accurate simulation !!!" 608 <<G4endl; 609 G4cout <<" *****************************************************************" 610 <<G4endl; 611 G4cout << G4endl; 612 613 return; 614 } 615 //////////////////////////////////////////////////////////////////////////////// 616 // 617