- Timestamp:
- Nov 5, 2010, 3:45:55 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/cascade/cascade/src/G4NucleiModel.cc
r1337 r1340 23 23 // * acceptance of all terms of the Geant4 Software license. * 24 24 // ******************************************************************** 25 // 26 // $Id: G4NucleiModel.cc,v 1.48.2.1 2010/06/25 09:44:56 gunter Exp $ 27 // Geant4 tag: $Name: geant4-09-04-beta-01 $ 25 // $Id: G4NucleiModel.cc,v 1.87 2010/09/24 06:26:06 mkelsey Exp $ 26 // Geant4 tag: $Name: geant4-09-03-ref-09 $ 28 27 // 29 28 // 20100112 M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly … … 41 40 // calculations. use G4Cascade*Channel for total xsec in pi-N 42 41 // N-N channels. Move absorptionCrossSection() from SpecialFunc. 43 44 //#define CHC_CHECK 42 // 20100610 M. Kelsey -- Replace another random-angle code block; add some 43 // diagnostic output for partner-list production. 44 // 20100617 M. Kelsey -- Replace preprocessor flag CHC_CHECK with 45 // G4CASCADE_DEBUG_CHARGE 46 // 20100620 M. Kelsey -- Improve error message on empty partners list, add 47 // four-momentum checking after EPCollider 48 // 20100621 M. Kelsey -- In boundaryTransition() account for momentum transfer 49 // to secondary by boosting into recoil nucleus "rest" frame. 50 // Don't need temporaries to create dummy "partners" for list. 51 // 20100622 M. Kelsey -- Restore use of "bindingEnergy()" function name, which 52 // is now a wrapper for G4NucleiProperties::GetBindingEnergy(). 53 // 20100623 M. Kelsey -- Eliminate some temporaries terminating partner-list, 54 // discard recoil boost for now. Initialize all data 55 // members in ctors. Allow generateModel() to be called 56 // mutliple times, by clearing vectors each time through; 57 // avoid extra work by returning if A and Z are same as 58 // before. 59 // 20100628 M. Kelsey -- Two momentum-recoil bugs; don't subtract energies! 60 // 20100715 M. Kelsey -- Make G4InuclNuclei in generateModel(), use for 61 // balance checking (only verbose>2) in generateParticleFate(). 62 // 20100721 M. Kelsey -- Use new G4CASCADE_CHECK_ECONS for conservation checks 63 // 20100723 M. Kelsey -- Move G4CollisionOutput buffer to .hh for reuse 64 // 20100726 M. Kelsey -- Preallocate arrays with number_of_zones dimension. 65 // 20100902 M. Kelsey -- Remove resize(3) directives from qdeutron/acsecs 66 // 20100907 M. Kelsey -- Limit interaction targets based on current nucleon 67 // counts (e.g., no pp if protonNumberCurrent < 2!) 68 // 20100914 M. Kelsey -- Migrate to integer A and Z 45 69 46 70 #include "G4NucleiModel.hh" 71 #include "G4CascadeCheckBalance.hh" 47 72 #include "G4CascadeInterpolator.hh" 48 73 #include "G4CascadeNNChannel.hh" … … 63 88 #include "G4LorentzConvertor.hh" 64 89 #include "G4Neutron.hh" 65 #include "G4NucleiProperties.hh"66 90 #include "G4Proton.hh" 67 91 … … 69 93 using namespace G4InuclSpecialFunctions; 70 94 71 72 95 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator; 73 96 74 G4NucleiModel::G4NucleiModel() : verboseLevel(0) { 75 if (verboseLevel > 3) { 76 G4cout << " >>> G4NucleiModel::G4NucleiModel" << G4endl; 77 } 78 } 79 80 G4NucleiModel::G4NucleiModel(G4InuclNuclei* nuclei) : verboseLevel(0) { 97 G4NucleiModel::G4NucleiModel() 98 : verboseLevel(0), A(0), Z(0), theNucleus(0), 99 neutronNumber(0), protonNumber(0), 100 neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0), 101 current_nucl2(0) {} 102 103 G4NucleiModel::G4NucleiModel(G4int a, G4int z) 104 : verboseLevel(0), A(0), Z(0), theNucleus(0), 105 neutronNumber(0), protonNumber(0), 106 neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0), 107 current_nucl2(0) { 108 generateModel(a,z); 109 } 110 111 G4NucleiModel::G4NucleiModel(G4InuclNuclei* nuclei) 112 : verboseLevel(0), A(0), Z(0), theNucleus(0), 113 neutronNumber(0), protonNumber(0), 114 neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0), 115 current_nucl2(0) { 116 generateModel(nuclei); 117 } 118 119 120 void 121 G4NucleiModel::generateModel(G4InuclNuclei* nuclei) { 81 122 generateModel(nuclei->getA(), nuclei->getZ()); 82 123 } … … 84 125 85 126 void 86 G4NucleiModel::generateModel(G4double a, G4double z) { 87 if (verboseLevel > 3) { 88 G4cout << " >>> G4NucleiModel::generateModel" << G4endl; 127 G4NucleiModel::generateModel(G4int a, G4int z) { 128 if (verboseLevel) { 129 G4cout << " >>> G4NucleiModel::generateModel A " << a << " Z " << z 130 << G4endl; 89 131 } 90 132 … … 95 137 const G4double pion_vp = 0.007; // in GeV 96 138 const G4double pion_vp_small = 0.007; 97 const G4double radForSmall = 8.0; // fermi98 const G4double piTimes4thirds = 4.189; // 4 Pi/3 139 const G4double radForSmall = 8.0; // Units 0.f fm? R_p = 0.8768 fm 140 const G4double piTimes4thirds = 4.189; // 4 Pi/3 FIXME! 99 141 const G4double mproton = G4Proton::Definition()->GetPDGMass() / GeV; 100 142 const G4double mneutron = G4Neutron::Definition()->GetPDGMass() / GeV; … … 102 144 const G4double alfa6[6] = { 0.9, 0.6, 0.4, 0.2, 0.1, 0.05 }; 103 145 146 // If model already built, just return; otherwise intialize everything 147 if (a == A && z == Z) { 148 if (verboseLevel > 1) 149 G4cout << " model already generated for A=" << a << ", Z=" << z << G4endl; 150 151 reset(); // Zeros out neutron/proton evaporates 152 return; 153 } 154 104 155 A = a; 105 156 Z = z; 157 delete theNucleus; 158 theNucleus = new G4InuclNuclei(A,Z); // For conservation checking 159 106 160 neutronNumber = a - z; 107 161 protonNumber = z; 108 neutronNumberCurrent = neutronNumber; 109 protonNumberCurrent = protonNumber; 110 111 // Set binding energies 112 // G4double dm = bindingEnergy(a, z); 113 G4double dm = G4NucleiProperties::GetBindingEnergy(G4lrint(a), G4lrint(z)); 114 115 binding_energies.push_back(0.001 * std::fabs(G4NucleiProperties::GetBindingEnergy(G4lrint(a-1), G4lrint(z-1)) - dm)); // for P 116 binding_energies.push_back(0.001 * std::fabs(G4NucleiProperties::GetBindingEnergy(G4lrint(a-1), G4lrint(z)) - dm)); // for N 162 reset(); 163 164 // Clear all parameters arrays for reloading 165 binding_energies.clear(); 166 nucleon_densities.clear(); 167 zone_potentials.clear(); 168 fermi_momenta.clear(); 169 zone_radii.clear(); 170 171 // Set binding energies 172 G4double dm = bindingEnergy(a,z); 173 174 binding_energies.push_back(0.001 * std::fabs(bindingEnergy(a-1,z-1)-dm)); // for P 175 binding_energies.push_back(0.001 * std::fabs(bindingEnergy(a-1,z)-dm)); // for N 117 176 118 177 G4double CU = cuu*G4cbrt(a); // half-density radius * 2.8197 … … 121 180 G4double CU2 = 0.0; 122 181 123 if (a > 4.5) { 124 std::vector<G4double> ur; 182 // This will be used to pre-allocate lots of arrays below 183 number_of_zones = (a < 5) ? 1 : (a < 100) ? 3 : 6; 184 185 // Buffers used for zone-by-zone calculations 186 G4double ur[7]; 187 G4double v[6]; 188 G4double v1[6]; 189 190 // These are passed into nested vectors, so can't be made C arrays 191 std::vector<G4double> rod(number_of_zones); 192 std::vector<G4double> pf(number_of_zones); 193 std::vector<G4double> vz(number_of_zones); 194 195 if (a > 4) { 125 196 G4int icase = 0; 126 197 127 if (a > 99.5) { 128 number_of_zones = 6; 129 ur.push_back(-D1); 130 198 if (a > 99) { 199 ur[0] = -D1; 131 200 for (G4int i = 0; i < number_of_zones; i++) { 132 201 G4double y = std::log((1.0 + D) / alfa6[i] - 1.0); 133 202 zone_radii.push_back(CU + AU * y); 134 ur .push_back(y);203 ur[i+1] = y; 135 204 } 136 137 } else if (a > 11.5) { 138 number_of_zones = 3; 139 ur.push_back(-D1); 140 205 } else if (a > 11) { 206 ur[0] = -D1; 141 207 for (G4int i = 0; i < number_of_zones; i++) { 142 208 G4double y = std::log((1.0 + D)/alfa3[i] - 1.0); 143 209 zone_radii.push_back(CU + AU * y); 144 ur .push_back(y);210 ur[i+1] = y; 145 211 } 146 147 212 } else { 148 number_of_zones = 3;149 213 icase = 1; 150 ur.push_back(0.0);151 214 152 215 G4double CU1 = CU * CU; 153 216 CU2 = std::sqrt(CU1 * (1.0 - 1.0 / a) + 6.4); 154 217 218 ur[0] = 0.0; 155 219 for (G4int i = 0; i < number_of_zones; i++) { 156 220 G4double y = std::sqrt(-std::log(alfa3[i])); 157 221 zone_radii.push_back(CU2 * y); 158 ur .push_back(y);222 ur[i+1] = y; 159 223 } 160 224 } 161 225 162 226 G4double tot_vol = 0.0; 163 std::vector<G4double> v;164 std::vector<G4double> v1;165 166 227 G4int i(0); 167 228 for (i = 0; i < number_of_zones; i++) { … … 169 230 170 231 if (icase == 0) { 171 v0 = volNumInt(ur[i], ur[i + 1], CU, D1); 172 232 v0 = volNumInt(ur[i], ur[i+1], CU, D1); 173 233 } else { 174 v0 = volNumInt1(ur[i], ur[i +1], CU2);175 } ;234 v0 = volNumInt1(ur[i], ur[i+1], CU2); 235 } 176 236 177 v .push_back(v0);237 v[i] = v0; 178 238 tot_vol += v0; 179 239 … … 181 241 if (i > 0) v0 -= zone_radii[i-1]*zone_radii[i-1]*zone_radii[i-1]; 182 242 183 v1 .push_back(v0);243 v1[i] = v0; 184 244 } 185 245 186 246 // Protons 187 G4double dd0 = z/tot_vol/piTimes4thirds;188 std::vector<G4double> rod;189 std::vector<G4double> pf;190 std::vector<G4double> vz;247 G4double dd0 = G4double(z)/tot_vol/piTimes4thirds; 248 rod.clear(); 249 pf.clear(); 250 vz.clear(); 191 251 192 252 for (i = 0; i < number_of_zones; i++) { … … 203 263 204 264 // Neutrons 205 dd0 = (a -z)/tot_vol/piTimes4thirds;265 dd0 = G4double(a-z)/tot_vol/piTimes4thirds; 206 266 rod.clear(); 207 267 pf.clear(); … … 214 274 pf.push_back(pff); 215 275 vz.push_back(0.5 * pff * pff / mneutron + binding_energies[1]); 216 } ;276 } 217 277 218 278 nucleon_densities.push_back(rod); … … 221 281 222 282 // pion stuff (primitive) 223 std::vector<G4double> vp(number_of_zones, pion_vp);283 const std::vector<G4double> vp(number_of_zones, pion_vp); 224 284 zone_potentials.push_back(vp); 225 285 226 286 // kaon potential (primitive) 227 std::vector<G4double> kp(number_of_zones, -0.015);287 const std::vector<G4double> kp(number_of_zones, -0.015); 228 288 zone_potentials.push_back(kp); 229 289 230 290 // hyperon potential (primitive) 231 std::vector<G4double> hp(number_of_zones, 0.03);291 const std::vector<G4double> hp(number_of_zones, 0.03); 232 292 zone_potentials.push_back(hp); 233 293 234 294 } else { // a < 5 235 number_of_zones = 1;236 295 G4double smallRad = radForSmall; 237 296 if (a == 4) smallRad *= 0.7; … … 240 299 241 300 // proton 242 std::vector<G4double> rod;243 std::vector<G4double> pf;244 std::vector<G4double> vz;245 301 for (G4int i = 0; i < number_of_zones; i++) { 246 302 G4double rd = vol*z; … … 272 328 273 329 // pion (primitive) 274 std::vector<G4double> vp(number_of_zones, pion_vp_small);330 const std::vector<G4double> vp(number_of_zones, pion_vp_small); 275 331 zone_potentials.push_back(vp); 276 332 277 333 // kaon potential (primitive) 278 std::vector<G4double> kp(number_of_zones, -0.015);334 const std::vector<G4double> kp(number_of_zones, -0.015); 279 335 zone_potentials.push_back(kp); 280 336 281 337 // hyperon potential (primitive) 282 std::vector<G4double> hp(number_of_zones, 0.03);338 const std::vector<G4double> hp(number_of_zones, 0.03); 283 339 zone_potentials.push_back(hp); 284 340 } 285 341 286 342 nuclei_radius = zone_radii[zone_radii.size() - 1]; 287 288 /*289 // Print nuclear radii and densities290 G4cout << " For A = " << a << " zone radii = ";291 for (G4int i = 0; i < number_of_zones; i++) G4cout << zone_radii[i] << " ";292 G4cout << " " << G4endl;293 294 G4cout << " proton densities: ";295 for (G4int i = 0; i < number_of_zones; i++)296 G4cout << nucleon_densities[0][i] << " ";297 G4cout << G4endl;298 299 G4cout << " neutron densities: ";300 for (G4int i = 0; i < number_of_zones; i++)301 G4cout << nucleon_densities[1][i] << " ";302 G4cout << G4endl;303 304 G4cout << " protons per shell " ;305 G4double rinner = 0.0;306 G4double router = 0.0;307 G4double shellVolume = 0.0;308 for (G4int i = 0; i < number_of_zones; i++) { // loop over zones309 router = zone_radii[i];310 shellVolume = piTimes4thirds*(router*router*router - rinner*rinner*rinner);311 G4cout << G4lrint(shellVolume*nucleon_densities[0][i]) << " ";312 rinner = router;313 }314 G4cout << G4endl;315 316 G4cout << " neutrons per shell " ;317 rinner = 0.0;318 router = 0.0;319 shellVolume = 0.0;320 for (G4int i = 0; i < number_of_zones; i++) { // loop over zones321 router = zone_radii[i];322 shellVolume = piTimes4thirds*(router*router*router - rinner*rinner*rinner);323 G4cout << G4lrint(shellVolume*nucleon_densities[1][i]) << " ";324 rinner = router;325 }326 G4cout << G4endl;327 */328 343 } 329 344 … … 345 360 G4NucleiModel::volNumInt(G4double r1, G4double r2, 346 361 G4double, G4double d1) const { 347 348 if (verboseLevel > 3) { 362 if (verboseLevel > 1) { 349 363 G4cout << " >>> G4NucleiModel::volNumInt" << G4endl; 350 364 } … … 373 387 r += dr1; 374 388 fi += r * (r + d2) / (1.0 + std::exp(r)); 375 } ;389 } 376 390 377 391 fun = 0.5 * fun1 + fi * dr; … … 398 412 G4NucleiModel::volNumInt1(G4double r1, G4double r2, 399 413 G4double cu2) const { 400 if (verboseLevel > 3) {414 if (verboseLevel > 1) { 401 415 G4cout << " >>> G4NucleiModel::volNumInt1" << G4endl; 402 416 } … … 447 461 448 462 void G4NucleiModel::printModel() const { 449 450 if (verboseLevel > 3) { 463 if (verboseLevel > 1) { 451 464 G4cout << " >>> G4NucleiModel::printModel" << G4endl; 452 465 } 453 466 454 467 G4cout << " nuclei model for A " << A << " Z " << Z << G4endl 455 << " proton binding energy " << binding_energies[0] <<456 457 << " N culei radius " << nuclei_radius << " number of zones " <<458 468 << " proton binding energy " << binding_energies[0] 469 << " neutron binding energy " << binding_energies[1] << G4endl 470 << " Nuclei radius " << nuclei_radius << " number of zones " 471 << number_of_zones << G4endl; 459 472 460 473 for (G4int i = 0; i < number_of_zones; i++) 461 462 474 G4cout << " zone " << i+1 << " radius " << zone_radii[i] << G4endl 463 475 << " protons: density " << getDensity(1,i) << " PF " << … … 480 492 G4InuclElementaryParticle 481 493 G4NucleiModel::generateNucleon(G4int type, G4int zone) const { 482 if (verboseLevel > 3) {494 if (verboseLevel > 1) { 483 495 G4cout << " >>> G4NucleiModel::generateNucleon" << G4endl; 484 496 } … … 492 504 G4NucleiModel::generateQuasiDeutron(G4int type1, G4int type2, 493 505 G4int zone) const { 494 495 if (verboseLevel > 3) { 506 if (verboseLevel > 1) { 496 507 G4cout << " >>> G4NucleiModel::generateQuasiDeutron" << G4endl; 497 508 } … … 518 529 void 519 530 G4NucleiModel::generateInteractionPartners(G4CascadParticle& cparticle) { 520 if (verboseLevel > 3) {531 if (verboseLevel > 1) { 521 532 G4cout << " >>> G4NucleiModel::generateInteractionPartners" << G4endl; 522 533 } 523 534 524 const G4double pi4by3 = 4.1887903; // 4 Pi / 3 535 const G4double pi4by3 = 4.1887903; // 4 Pi / 3 FIXME! 525 536 const G4double small = 1.0e-10; 526 537 const G4double huge_num = 50.0; … … 556 567 r_in = zone_radii[zone - 1]; 557 568 r_out = zone_radii[zone]; 558 } ;569 } 559 570 560 571 G4double path = cparticle.getPathToTheNextZone(r_in, r_out); 561 572 562 if (verboseLevel > 2) {573 if (verboseLevel > 2) { 563 574 G4cout << " r_in " << r_in << " r_out " << r_out << " path " << path << G4endl; 564 575 } 565 576 566 if (path < -small) { // something wrong 577 if (path < -small) { // something wrong 578 G4cerr << " generateInteractionPartners-> negative path length" << G4endl; 567 579 return; 568 569 } else if (std::fabs(path) < small) { // just on the boundary 570 path = 0.0; 571 572 G4InuclElementaryParticle particle; // Dummy -- no type or momentum 573 thePartners.push_back(partner(particle, path)); 574 575 } else { // normal case 576 G4LorentzConvertor dummy_convertor; 577 dummy_convertor.setBullet(pmom, pmass); 578 579 for (G4int ip = 1; ip < 3; ip++) { 580 G4InuclElementaryParticle particle = generateNucleon(ip, zone); 581 dummy_convertor.setTarget(particle.getMomentum(), particle.getMass()); 580 } 581 582 if (std::fabs(path) < small) { // just on the boundary 583 if (verboseLevel > 3) 584 G4cout << " generateInteractionPartners-> zero path" << G4endl; 585 586 thePartners.push_back(partner()); // Dummy list terminator with zero path 587 return; 588 } 589 590 G4LorentzConvertor dummy_convertor; 591 dummy_convertor.setBullet(pmom, pmass); 592 593 for (G4int ip = 1; ip < 3; ip++) { 594 // Only process nucleons which remain active in target 595 if (ip==1 && protonNumberCurrent < 1) continue; 596 if (ip==2 && neutronNumberCurrent < 1) continue; 597 598 // All nucleons are assumed to be at rest when colliding 599 G4InuclElementaryParticle particle = generateNucleon(ip, zone); 600 dummy_convertor.setTarget(particle.getMomentum(), particle.getMass()); 601 G4double ekin = dummy_convertor.getKinEnergyInTheTRS(); 602 603 G4double csec = totalCrossSection(ekin, ptype * ip); 604 605 if(verboseLevel > 2) { 606 G4cout << " ip " << ip << " ekin " << ekin << " csec " << csec << G4endl; 607 } 608 609 G4double dens = nucleon_densities[ip - 1][zone]; 610 G4double rat = getRatio(ip); 611 G4double pw = -path * dens * csec * rat; 612 613 if (pw < -huge_num) pw = -huge_num; 614 pw = 1.0 - std::exp(pw); 615 616 if (verboseLevel > 2) { 617 G4cout << " pw " << pw << " rat " << rat << G4endl; 618 } 619 620 G4double spath = path; 621 622 if (inuclRndm() < pw) { 623 spath = -1.0 / dens / csec / rat * std::log(1.0 - pw * inuclRndm()); 624 if (cparticle.young(young_cut, spath)) spath = path; 625 626 if (verboseLevel > 2) { 627 G4cout << " ip " << ip << " spath " << spath << G4endl; 628 } 629 } 630 631 if (spath < path) { 632 if (verboseLevel > 3) { 633 G4cout << " adding partner[" << thePartners.size() << "]: "; 634 particle.printParticle(); 635 } 636 thePartners.push_back(partner(particle, spath)); 637 } 638 } 639 640 if (verboseLevel > 2) { 641 G4cout << " after nucleons " << thePartners.size() << " path " << path << G4endl; 642 } 643 644 if (cparticle.getParticle().pion()) { // absorption possible 645 if (verboseLevel > 2) { 646 G4cout << " trying quasi-deuterons with bullet: "; 647 cparticle.getParticle().printParticle(); 648 } 649 650 // Initialize buffers for quasi-deuteron results 651 qdeutrons.clear(); 652 acsecs.clear(); 653 654 G4double tot_abs_csec = 0.0; 655 G4double abs_sec; 656 G4double vol = zone_radii[zone]*zone_radii[zone]*zone_radii[zone]; 657 658 if (zone > 0) vol -= zone_radii[zone-1]*zone_radii[zone-1]*zone_radii[zone-1]; 659 vol *= pi4by3; 660 661 G4double rat = getRatio(1); 662 G4double rat1 = getRatio(2); 663 664 // FIXME: Shouldn't be creating zero-cross-section states! 665 666 // Proton-proton state interacts with pi-, pi0 only 667 G4InuclElementaryParticle ppd = generateQuasiDeutron(1, 1, zone); 668 669 if (protonNumberCurrent < 2 || !(ptype == pi0 || ptype == pim)) { 670 abs_sec = 0.0; 671 } else { 672 dummy_convertor.setTarget(ppd.getMomentum(), ppd.getMass()); 673 582 674 G4double ekin = dummy_convertor.getKinEnergyInTheTRS(); 583 584 // Total cross section converted from mb to fm**2 585 G4double csec = totalCrossSection(ekin, ptype * ip); 586 587 if(verboseLevel > 2){ 588 G4cout << " ip " << ip << " ekin " << ekin << " csec " << csec << G4endl; 675 676 if (verboseLevel > 2) { 677 G4cout << " ptype=" << ptype << " using pp target" << G4endl; 678 ppd.printParticle(); 589 679 } 590 591 G4double dens = nucleon_densities[ip - 1][zone]; 592 G4double rat = getRatio(ip); 593 G4double pw = -path * dens * csec * rat; 594 595 if (pw < -huge_num) pw = -huge_num; 596 pw = 1.0 - std::exp(pw); 597 598 if (verboseLevel > 2){ 599 G4cout << " pw " << pw << " rat " << rat << G4endl; 600 } 601 602 G4double spath = path; 603 604 if (inuclRndm() < pw) { 605 spath = -1.0 / dens / csec / rat * std::log(1.0 - pw * inuclRndm()); 606 if (cparticle.young(young_cut, spath)) spath = path; 607 608 if (verboseLevel > 2){ 609 G4cout << " ip " << ip << " spath " << spath << G4endl; 610 } 611 612 }; 613 if (spath < path) thePartners.push_back(partner(particle, spath)); 614 }; 615 616 if (verboseLevel > 2){ 617 G4cout << " after nucleons " << thePartners.size() << " path " << path << G4endl; 618 } 619 620 if (cparticle.getParticle().pion()) { // absorption possible 621 if (verboseLevel > 2) { 622 G4cout << " trying quasi-deuterons with bullet: "; 623 cparticle.getParticle().printParticle(); 624 } 625 626 std::vector<G4InuclElementaryParticle> qdeutrons(3); 627 std::vector<G4double> acsecs(3); 628 629 G4double tot_abs_csec = 0.0; 630 G4double abs_sec; 631 G4double vol = zone_radii[zone]*zone_radii[zone]*zone_radii[zone]; 632 633 if (zone > 0) vol -= zone_radii[zone-1]*zone_radii[zone-1]*zone_radii[zone-1]; 634 vol *= pi4by3; 635 636 G4double rat = getRatio(1); 637 G4double rat1 = getRatio(2); 638 639 G4InuclElementaryParticle ppd = generateQuasiDeutron(1, 1, zone); 640 641 if (ptype == 7 || ptype == 5) { 642 dummy_convertor.setTarget(ppd.getMomentum(), ppd.getMass()); 643 644 G4double ekin = dummy_convertor.getKinEnergyInTheTRS(); 645 646 if (verboseLevel > 2) { 647 G4cout << " ptype=" << ptype << " using pp target" << G4endl; 648 ppd.printParticle(); 649 } 650 651 abs_sec = absorptionCrossSection(ekin, ptype); 652 abs_sec *= nucleon_densities[0][zone] * nucleon_densities[0][zone]* 653 rat * rat * vol; 654 655 } else { 656 abs_sec = 0.0; 657 }; 658 659 // abs_sec = 0.0; 660 tot_abs_csec += abs_sec; 661 acsecs.push_back(abs_sec); 662 qdeutrons.push_back(ppd); 663 664 G4InuclElementaryParticle npd = generateQuasiDeutron(1, 2, zone); 665 680 681 abs_sec = absorptionCrossSection(ekin, ptype); 682 abs_sec *= nucleon_densities[0][zone] * nucleon_densities[0][zone]* 683 rat * rat * vol; 684 } 685 686 tot_abs_csec += abs_sec; 687 acsecs.push_back(abs_sec); 688 qdeutrons.push_back(ppd); 689 690 // Proton-neutron state interacts with any pion type 691 G4InuclElementaryParticle npd = generateQuasiDeutron(1, 2, zone); 692 693 if (protonNumberCurrent < 1 || neutronNumberCurrent < 1) { 694 abs_sec = 0.0; 695 } else { 666 696 dummy_convertor.setTarget(npd.getMomentum(), npd.getMass()); 667 697 668 698 G4double ekin = dummy_convertor.getKinEnergyInTheTRS(); 669 699 670 700 if (verboseLevel > 2) { 671 701 G4cout << " using np target" << G4endl; 672 702 npd.printParticle(); 673 703 } 674 704 675 705 abs_sec = absorptionCrossSection(ekin, ptype); 676 706 abs_sec *= pn_spec * nucleon_densities[0][zone] * nucleon_densities[1][zone] * 677 rat * rat1 * vol; 678 tot_abs_csec += abs_sec; 679 acsecs.push_back(abs_sec); 680 qdeutrons.push_back(npd); 681 682 G4InuclElementaryParticle nnd = generateQuasiDeutron(2, 2, zone); 683 684 if (ptype == 7 || ptype == 3) { 685 dummy_convertor.setTarget(nnd.getMomentum(), nnd.getMass()); 686 687 G4double ekin = dummy_convertor.getKinEnergyInTheTRS(); 688 689 if (verboseLevel > 2) { 690 G4cout << " ptype=" << ptype << " using nn target" << G4endl; 691 nnd.printParticle(); 707 rat * rat1 * vol; 708 } 709 710 tot_abs_csec += abs_sec; 711 acsecs.push_back(abs_sec); 712 qdeutrons.push_back(npd); 713 714 // Neutron-neutron state interacts with pi+, pi0 only 715 G4InuclElementaryParticle nnd = generateQuasiDeutron(2, 2, zone); 716 717 if (neutronNumberCurrent < 2 || !(ptype == pi0 || ptype == pip)) { 718 abs_sec = 0.0; 719 } else { 720 dummy_convertor.setTarget(nnd.getMomentum(), nnd.getMass()); 721 722 G4double ekin = dummy_convertor.getKinEnergyInTheTRS(); 723 724 if (verboseLevel > 2) { 725 G4cout << " ptype=" << ptype << " using nn target" << G4endl; 726 nnd.printParticle(); 727 } 728 729 abs_sec = absorptionCrossSection(ekin, ptype); 730 abs_sec *= nucleon_densities[1][zone] * nucleon_densities[1][zone] * 731 rat1 * rat1 * vol; 732 } 733 734 tot_abs_csec += abs_sec; 735 acsecs.push_back(abs_sec); 736 qdeutrons.push_back(nnd); 737 738 // Select quasideuteron interaction from non-zero cross-section choices 739 if (verboseLevel > 2){ 740 G4cout << " rod1 " << acsecs[0] << " rod2 " << acsecs[1] 741 << " rod3 " << acsecs[2] << G4endl; 742 } 743 744 if (tot_abs_csec > small) { 745 G4double pw = -path * tot_abs_csec; 746 747 if (pw < -huge_num) pw = -huge_num; 748 pw = 1.0 - std::exp(pw); 749 750 if (verboseLevel > 2){ 751 G4cout << " pw " << pw << G4endl; 752 } 753 754 G4double apath = path; 755 756 if (inuclRndm() < pw) 757 apath = -1.0 / tot_abs_csec * std::log(1.0 - pw * inuclRndm()); 758 759 if (cparticle.young(young_cut, apath)) apath = path; 760 761 if(verboseLevel > 2){ 762 G4cout << " apath " << apath << " path " << path << G4endl; 763 } 764 765 if (apath < path) { // choose the qdeutron 766 G4double sl = inuclRndm() * tot_abs_csec; 767 G4double as = 0.0; 768 769 for (G4int i = 0; i < 3; i++) { 770 as += acsecs[i]; 771 if (sl < as) { 772 if (verboseLevel > 2) G4cout << " deut type " << i << G4endl; 773 thePartners.push_back(partner(qdeutrons[i], apath)); 774 break; 775 } 692 776 } 693 694 abs_sec = absorptionCrossSection(ekin, ptype); 695 abs_sec *= nucleon_densities[1][zone] * nucleon_densities[1][zone] * 696 rat1 * rat1 * vol; 697 698 } else { 699 abs_sec = 0.0; 700 }; 701 702 // abs_sec = 0.0; 703 tot_abs_csec += abs_sec; 704 acsecs.push_back(abs_sec); 705 qdeutrons.push_back(nnd); 706 707 if (verboseLevel > 2){ 708 G4cout << " rod1 " << acsecs[0] << " rod2 " << acsecs[1] 709 << " rod3 " << acsecs[2] << G4endl; 710 } 711 712 if (tot_abs_csec > small) { 713 714 G4double pw = -path * tot_abs_csec; 715 716 if (pw < -huge_num) pw = -huge_num; 717 pw = 1.0 - std::exp(pw); 718 719 if (verboseLevel > 2){ 720 G4cout << " pw " << pw << G4endl; 721 } 722 723 G4double apath = path; 724 725 if (inuclRndm() < pw) 726 apath = -1.0 / tot_abs_csec * std::log(1.0 - pw * inuclRndm()); 727 728 if (cparticle.young(young_cut, apath)) apath = path; 729 730 if(verboseLevel > 2){ 731 G4cout << " apath " << apath << " path " << path << G4endl; 732 } 733 734 if (apath < path) { // chose the qdeutron 735 736 G4double sl = inuclRndm() * tot_abs_csec; 737 G4double as = 0.0; 738 739 for (G4int i = 0; i < 3; i++) { 740 as += acsecs[i]; 741 742 if (sl < as) { 743 744 if (verboseLevel > 2){ 745 G4cout << " deut type " << i << G4endl; 746 } 747 748 thePartners.push_back(partner(qdeutrons[i], apath)); 749 750 break; 751 }; 752 }; 753 }; 754 }; 755 }; 756 757 if(verboseLevel > 2){ 758 G4cout << " after deutrons " << thePartners.size() << G4endl; 759 } 760 761 if (thePartners.size() > 1) { // Sort list by path length 762 std::sort(thePartners.begin(), thePartners.end(), sortPartners); 763 } 764 765 G4InuclElementaryParticle particle; // Dummy for end of list 766 thePartners.push_back(partner(particle, path)); 767 } 768 769 return; 777 } 778 } 779 } 780 781 if (verboseLevel > 2) { 782 G4cout << " after deuterons " << thePartners.size() << " partners" 783 << G4endl; 784 } 785 786 if (thePartners.size() > 1) { // Sort list by path length 787 std::sort(thePartners.begin(), thePartners.end(), sortPartners); 788 } 789 790 G4InuclElementaryParticle particle; // Total path at end of list 791 thePartners.push_back(partner(particle, path)); 770 792 } 771 793 … … 774 796 G4NucleiModel::generateParticleFate(G4CascadParticle& cparticle, 775 797 G4ElementaryParticleCollider* theElementaryParticleCollider) { 776 if (verboseLevel > 3)798 if (verboseLevel > 1) 777 799 G4cout << " >>> G4NucleiModel::generateParticleFate" << G4endl; 778 800 801 if (verboseLevel > 2) { 802 G4cout << " cparticle: "; 803 cparticle.print(); 804 } 805 806 // Create four-vector checking 807 #ifdef G4CASCADE_CHECK_ECONS 808 G4CascadeCheckBalance balance(0.005, 0.01, "G4NucleiModel"); // Second arg is in GeV 809 balance.setVerboseLevel(verboseLevel); 810 #endif 811 779 812 outgoing_cparticles.clear(); // Clear return buffer for this event 780 781 813 generateInteractionPartners(cparticle); // Fills "thePartners" data 782 814 783 815 if (thePartners.empty()) { // smth. is wrong -> needs special treatment 784 G4cout << " generateParticleFate-> can not be here " << G4endl; 816 G4cerr << " generateParticleFate-> got empty interaction-partners list " 817 << G4endl; 785 818 return outgoing_cparticles; 786 819 } 787 820 788 G4int npart = thePartners.size(); 789 790 if (npart == 1) { // cparticle is on the next zone entry 791 // need to go here if particle outside nucleus ? 792 // 821 G4int npart = thePartners.size(); // Last item is a total-path placeholder 822 823 if (npart == 1) { // cparticle is on the next zone entry 793 824 cparticle.propagateAlongThePath(thePartners[0].second); 794 825 cparticle.incrementCurrentPath(thePartners[0].second); … … 796 827 outgoing_cparticles.push_back(cparticle); 797 828 798 if (verboseLevel > 2) {829 if (verboseLevel > 2) { 799 830 G4cout << " next zone " << G4endl; 800 831 cparticle.print(); 801 832 } 802 803 } else { // there are possible interactions 804 833 } else { // there are possible interactions 834 if (verboseLevel > 1) 835 G4cout << " processing " << npart-1 << " possible interactions" << G4endl; 836 805 837 G4ThreeVector old_position = cparticle.getPosition(); 806 807 G4InuclElementaryParticle bullet = cparticle.getParticle(); 808 838 G4InuclElementaryParticle& bullet = cparticle.getParticle(); 809 839 G4bool no_interaction = true; 810 811 840 G4int zone = cparticle.getCurrentZone(); 812 841 813 G4CollisionOutput output; 814 815 for (G4int i = 0; i < npart - 1; i++) { 842 for (G4int i=0; i<npart-1; i++) { // Last item is a total-path placeholder 816 843 if (i > 0) cparticle.updatePosition(old_position); 817 844 818 G4InuclElementaryParticle target = thePartners[i].first;819 820 if (verboseLevel > 2){821 if (target.quasi_deutron()) 822 G4cout << " try absorption: target " << target.type() << " bullet " <<823 bullet.type()<< G4endl;845 G4InuclElementaryParticle& target = thePartners[i].first; 846 847 if (verboseLevel > 3) { 848 if (target.quasi_deutron()) G4cout << " try absorption: "; 849 G4cout << " target " << target.type() << " bullet " << bullet.type() 850 << G4endl; 824 851 } 825 852 826 output.reset(); 827 theElementaryParticleCollider->collide(&bullet, &target, output); 828 829 if (verboseLevel > 2) output.printCollisionOutput(); 853 EPCoutput.reset(); 854 theElementaryParticleCollider->collide(&bullet, &target, EPCoutput); 855 856 if (verboseLevel > 2) { 857 EPCoutput.printCollisionOutput(); 858 #ifdef G4CASCADE_CHECK_ECONS 859 balance.collide(&bullet, &target, EPCoutput); 860 balance.okay(); // Do checks, but ignore result 861 #endif 862 } 830 863 831 864 // Don't need to copy list, as "output" isn't changed again below 832 865 const std::vector<G4InuclElementaryParticle>& outgoing_particles = 833 output.getOutgoingParticles(); 834 835 if (passFermi(outgoing_particles, zone)) { // interaction 836 cparticle.propagateAlongThePath(thePartners[i].second); 837 G4ThreeVector new_position = cparticle.getPosition(); 838 839 for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++) { 840 G4CascadParticle temp(outgoing_particles[ip], new_position, zone, 0.0, 0); 841 outgoing_cparticles.push_back(temp); 842 } 843 844 no_interaction = false; 845 current_nucl1 = 0; 846 current_nucl2 = 0; 847 #ifdef CHC_CHECK 866 EPCoutput.getOutgoingParticles(); 867 868 if (!passFermi(outgoing_particles, zone)) continue; // Interaction fails 869 870 // Successful interaction, add results to output list 871 cparticle.propagateAlongThePath(thePartners[i].second); 872 G4ThreeVector new_position = cparticle.getPosition(); 873 874 if (verboseLevel > 2) 875 G4cout << " adding " << outgoing_particles.size() 876 << " output particles" << G4endl; 877 878 for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++) { 879 G4CascadParticle temp(outgoing_particles[ip], new_position, zone, 0.0, 0); 880 outgoing_cparticles.push_back(temp); 881 } 882 883 no_interaction = false; 884 current_nucl1 = 0; 885 current_nucl2 = 0; 886 887 #ifdef G4CASCADE_DEBUG_CHARGE 888 { 848 889 G4double out_charge = 0.0; 849 890 … … 851 892 out_charge += outgoing_particles[ip].getCharge(); 852 893 853 G4cout << " multiplicity " << outgoing_particles.size() << 854 " bul type " << bullet.type() << " targ type " << target.type() << 855 G4endl << " initial charge " << bullet.getCharge() + target.getCharge() 894 G4cout << " multiplicity " << outgoing_particles.size() 895 << " bul type " << bullet.type() 896 << " targ type " << target.type() 897 << "\n initial charge " 898 << bullet.getCharge() + target.getCharge() 856 899 << " out charge " << out_charge << G4endl; 900 } 857 901 #endif 858 859 if (verboseLevel > 2){ 860 G4cout << " partner type " << target.type() << G4endl; 861 } 862 863 if (target.nucleon()) { 864 current_nucl1 = target.type(); 865 866 } else { 867 if (verboseLevel > 2) G4cout << " good absorption " << G4endl; 868 869 current_nucl1 = (target.type() - 100) / 10; 870 current_nucl2 = target.type() - 100 - 10 * current_nucl1; 871 } 872 873 if (current_nucl1 == 1) { 874 protonNumberCurrent -= 1.0; 875 876 } else { 877 neutronNumberCurrent -= 1.0; 878 }; 879 880 if (current_nucl2 == 1) { 881 protonNumberCurrent -= 1.0; 882 883 } else if(current_nucl2 == 2) { 884 neutronNumberCurrent -= 1.0; 885 }; 886 887 break; 888 }; 889 } // loop over partners 890 891 if (no_interaction) { // still no interactions 902 903 if (verboseLevel > 2) 904 G4cout << " partner type " << target.type() << G4endl; 905 906 if (target.nucleon()) { 907 current_nucl1 = target.type(); 908 } else { 909 if (verboseLevel > 2) G4cout << " good absorption " << G4endl; 910 current_nucl1 = (target.type() - 100) / 10; 911 current_nucl2 = target.type() - 100 - 10 * current_nucl1; 912 } 913 914 if (current_nucl1 == 1) { 915 if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl; 916 protonNumberCurrent--; 917 } else { 918 if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl; 919 neutronNumberCurrent--; 920 } 921 922 if (current_nucl2 == 1) { 923 if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl; 924 protonNumberCurrent--; 925 } else if (current_nucl2 == 2) { 926 if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl; 927 neutronNumberCurrent--; 928 } 929 930 break; 931 } // loop over partners 932 933 if (no_interaction) { // still no interactions 934 if (verboseLevel > 1) G4cout << " no interaction " << G4endl; 935 936 // For conservation checking (below), get particle before updating 937 static G4InuclElementaryParticle prescatCP; // Avoid memory churn 938 prescatCP = cparticle.getParticle(); 939 940 // Last "partner" is just a total-path placeholder 892 941 cparticle.updatePosition(old_position); 893 cparticle.propagateAlongThePath(thePartners[npart -1].second);894 cparticle.incrementCurrentPath(thePartners[npart -1].second);942 cparticle.propagateAlongThePath(thePartners[npart-1].second); 943 cparticle.incrementCurrentPath(thePartners[npart-1].second); 895 944 boundaryTransition(cparticle); 896 945 outgoing_cparticles.push_back(cparticle); 897 }; 898 }; 946 947 // Check conservation for simple scattering (ignore target nucleus!) 948 #ifdef G4CASCADE_CHECK_ECONS 949 if (verboseLevel > 2) { 950 balance.collide(&prescatCP, 0, outgoing_cparticles); 951 balance.okay(); // Report violations, but don't act on them 952 } 953 #endif 954 } 955 } // if (npart == 1) [else] 899 956 900 957 return outgoing_cparticles; … … 903 960 G4bool G4NucleiModel::passFermi(const std::vector<G4InuclElementaryParticle>& particles, 904 961 G4int zone) { 962 if (verboseLevel > 1) { 963 G4cout << " >>> G4NucleiModel::passFermi" << G4endl; 964 } 965 966 // Only check Fermi momenta for nucleons 967 for (G4int i = 0; i < G4int(particles.size()); i++) { 968 if (!particles[i].nucleon()) continue; 969 970 G4int type = particles[i].type(); 971 G4double mom = particles[i].getMomModule(); 972 G4double pf = fermi_momenta[type-1][zone]; 973 974 if (verboseLevel > 2) 975 G4cout << " type " << type << " p " << mom << " pf " << pf << G4endl; 976 977 if (mom < pf) { 978 if (verboseLevel > 2) G4cout << " rejected by Fermi" << G4endl; 979 return false; 980 } 981 } 982 return true; 983 } 984 985 void G4NucleiModel::boundaryTransition(G4CascadParticle& cparticle) { 986 if (verboseLevel > 1) { 987 G4cout << " >>> G4NucleiModel::boundaryTransition" << G4endl; 988 } 989 990 G4int zone = cparticle.getCurrentZone(); 991 992 if (cparticle.movingInsideNuclei() && zone == 0) { 993 G4cerr << " boundaryTransition-> in zone 0 " << G4endl; 994 return; 995 } 996 997 G4LorentzVector mom = cparticle.getMomentum(); 998 G4ThreeVector pos = cparticle.getPosition(); 999 1000 G4int type = cparticle.getParticle().type(); 1001 1002 G4double pr = pos.dot(mom.vect()); 1003 G4double r = pos.mag(); 1004 1005 pr /= r; 1006 1007 G4int next_zone = cparticle.movingInsideNuclei() ? zone - 1 : zone + 1; 1008 1009 G4double dv = getPotential(type,zone) - getPotential(type, next_zone); 1010 // G4cout << "Potentials for type " << type << " = " 1011 // << getPotential(type,zone) << " , " 1012 // << getPotential(type,next_zone) << G4endl; 1013 1014 G4double qv = dv * dv - 2.0 * dv * mom.e() + pr * pr; 1015 1016 G4double p1r; 1017 905 1018 if (verboseLevel > 3) { 906 G4cout << " >>> G4NucleiModel::passFermi" << G4endl; 907 } 908 909 for (G4int i = 0; i < G4int(particles.size()); i++) { 910 911 if (particles[i].nucleon()) { 912 913 if (verboseLevel > 2){ 914 G4cout << " type " << particles[i].type() << " p " << particles[i].getMomModule() 915 << " pf " << fermi_momenta[particles[i].type() - 1][zone] << G4endl; 916 } 917 918 if (particles[i].getMomModule() < fermi_momenta[particles[i].type() - 1][zone]) { 919 920 if (verboseLevel > 2) { 921 G4cout << " rejected by fermi: type " << particles[i].type() << 922 " p " << particles[i].getMomModule() << G4endl; 923 } 924 925 return false; 926 }; 927 }; 928 }; 929 return true; 930 } 931 932 void G4NucleiModel::boundaryTransition(G4CascadParticle& cparticle) { 933 1019 G4cout << " type " << type << " zone " << zone << " next " << next_zone 1020 << " qv " << qv << " dv " << dv << G4endl; 1021 } 1022 1023 if (qv <= 0.0) { // reflection 1024 if (verboseLevel > 3) G4cout << " reflects off boundary" << G4endl; 1025 p1r = -pr; 1026 cparticle.incrementReflectionCounter(); 1027 } else { // transition 1028 if (verboseLevel > 3) G4cout << " passes thru boundary" << G4endl; 1029 p1r = std::sqrt(qv); 1030 if(pr < 0.0) p1r = -p1r; 1031 cparticle.updateZone(next_zone); 1032 cparticle.resetReflection(); 1033 } 1034 1035 G4double prr = (p1r - pr) / r; 1036 934 1037 if (verboseLevel > 3) { 935 G4cout << " >>> G4NucleiModel::boundaryTransition" << G4endl; 936 } 937 938 G4int zone = cparticle.getCurrentZone(); 939 940 if (cparticle.movingInsideNuclei() && zone == 0) { 941 G4cout << " boundaryTransition-> in zone 0 " << G4endl; 942 943 } else { 944 G4LorentzVector mom = cparticle.getMomentum(); 945 G4ThreeVector pos = cparticle.getPosition(); 946 947 G4int type = cparticle.getParticle().type(); 948 949 G4double pr = pos.dot(mom.vect()); 950 G4double r = pos.mag(); 951 952 pr /= r; 953 954 G4int next_zone = cparticle.movingInsideNuclei() ? zone - 1 : zone + 1; 955 956 G4double dv = getPotential(type,zone) - getPotential(type, next_zone); 957 // G4cout << "Potentials for type " << type << " = " 958 // << getPotential(type,zone) << " , " 959 // << getPotential(type,next_zone) << G4endl; 960 961 G4double qv = dv * dv - 2.0 * dv * mom.e() + pr * pr; 962 963 G4double p1r; 964 965 if (verboseLevel > 2){ 966 G4cout << " type " << type << " zone " << zone 967 << " next " << next_zone 968 << " qv " << qv << " dv " << dv << G4endl; 969 } 970 971 if(qv <= 0.0) { // reflection 972 p1r = -pr; 973 cparticle.incrementReflectionCounter(); 974 975 } else { // transition 976 p1r = std::sqrt(qv); 977 if(pr < 0.0) p1r = -p1r; 978 cparticle.updateZone(next_zone); 979 cparticle.resetReflection(); 980 }; 981 982 G4double prr = (p1r - pr) / r; 983 984 mom.setVect(mom.vect() + pos*prr); 985 cparticle.updateParticleMomentum(mom); 986 }; 1038 G4cout << " prr " << prr << " delta px " << prr*pos.x() << " py " 1039 << prr*pos.y() << " pz " << prr*pos.z() << " mag " 1040 << std::fabs(prr*r) << G4endl; 1041 } 1042 1043 mom.setVect(mom.vect() + pos*prr); 1044 cparticle.updateParticleMomentum(mom); 987 1045 } 988 1046 989 1047 G4bool G4NucleiModel::worthToPropagate(const G4CascadParticle& cparticle) const { 990 991 if (verboseLevel > 3) { 1048 if (verboseLevel > 1) { 992 1049 G4cout << " >>> G4NucleiModel::worthToPropagate" << G4endl; 993 1050 } … … 1010 1067 << worth << G4endl; 1011 1068 } 1012 } ;1069 } 1013 1070 1014 1071 return worth; 1015 1072 } 1016 1073 1074 1017 1075 G4double G4NucleiModel::getRatio(G4int ip) const { 1018 1019 if (verboseLevel > 3) { 1076 if (verboseLevel > 1) { 1020 1077 G4cout << " >>> G4NucleiModel::getRatio" << G4endl; 1021 1078 } 1022 1079 1023 G4double rat;1024 // G4double ratm;1025 1026 // Calculate number of protons and neutrons in local region1027 // G4double Athird = G4cbrt(A);1028 // G4double Nneut = Athird*(A-Z)/A;1029 // G4double Nprot = Athird*Z/A;1030 1031 // Reduce number of1032 1080 if (ip == 1) { 1081 if (verboseLevel > 2) { 1082 G4cout << " current " << protonNumberCurrent << " inp " << protonNumber 1083 << G4endl; 1084 } 1085 return G4double(protonNumberCurrent)/G4double(protonNumber); 1086 } 1087 1088 if (ip == 2) { 1033 1089 if (verboseLevel > 2){ 1034 G4cout << " current " << protonNumberCurrent << " inp " << protonNumber << G4endl; 1035 } 1036 1037 rat = protonNumberCurrent/protonNumber; 1038 1039 // Calculate ratio modified for local region 1040 // G4double deltaP = protonNumber - protonNumberCurrent; 1041 // G4cout << " deltaP = " << deltaP << G4endl; 1042 // ratm = std::max(0.0, (Nprot - deltaP)/Nprot); 1043 1044 } else { 1045 if (verboseLevel > 2){ 1046 G4cout << " current " << neutronNumberCurrent << " inp " << neutronNumber << G4endl; 1047 } 1048 1049 rat = neutronNumberCurrent/neutronNumber; 1050 1051 // Calculate ratio modified for local region 1052 // G4double deltaN = neutronNumber - neutronNumberCurrent; 1053 // G4cout << " deltaN = " << deltaN << G4endl; 1054 // ratm = std::max(0.0, (Nneut - deltaN)/Nneut); 1055 } 1056 1057 // G4cout << " get ratio: ratm = " << ratm << G4endl; 1058 return rat; 1059 // return ratm; 1060 } 1061 1062 G4CascadParticle G4NucleiModel::initializeCascad(G4InuclElementaryParticle* particle) { 1063 1064 if (verboseLevel > 3) { 1065 G4cout << " >>> G4NucleiModel::initializeCascad(G4InuclElementaryParticle* particle)" << G4endl; 1090 G4cout << " current " << neutronNumberCurrent << " inp " << neutronNumber 1091 << G4endl; 1092 } 1093 return G4double(neutronNumberCurrent)/G4double(neutronNumber); 1094 } 1095 1096 return 0.; 1097 } 1098 1099 G4CascadParticle 1100 G4NucleiModel::initializeCascad(G4InuclElementaryParticle* particle) { 1101 if (verboseLevel > 1) { 1102 G4cout << " >>> G4NucleiModel::initializeCascad(particle)" << G4endl; 1066 1103 } 1067 1104 1068 1105 const G4double large = 1000.0; 1069 1106 1070 G4double s1 = std::sqrt(inuclRndm()); 1071 G4double phi = randomPHI(); 1072 G4double rz = nuclei_radius * s1; 1073 1074 G4ThreeVector pos(rz*std::cos(phi), rz*std::sin(phi), 1075 -nuclei_radius*std::sqrt(1.0 - s1*s1)); 1076 1107 // FIXME: Previous version generated random sin(theta), then used -cos(theta) 1108 // Using generateWithRandomAngles changes result! 1109 // G4ThreeVector pos = generateWithRandomAngles(nuclei_radius).vect(); 1110 G4double costh = std::sqrt(1.0 - inuclRndm()); 1111 G4ThreeVector pos = generateWithFixedTheta(-costh, nuclei_radius); 1112 1077 1113 G4CascadParticle cpart(*particle, pos, number_of_zones, large, 0); 1078 1114 … … 1085 1121 G4InuclNuclei* target, 1086 1122 modelLists& output) { 1087 1088 if (verboseLevel > 3) {1089 G4cout << " >>> G4NucleiModel::initializeCascad(G4InuclNuclei* bullet, G4InuclNuclei* target)"<< G4endl;1123 if (verboseLevel) { 1124 G4cout << " >>> G4NucleiModel::initializeCascad(bullet,target,output)" 1125 << G4endl; 1090 1126 } 1091 1127 … … 1110 1146 1111 1147 // first decide whether it will be cascad or compound final nuclei 1112 G4 doubleab = bullet->getA();1113 G4 doublezb = bullet->getZ();1114 G4 doubleat = target->getA();1115 G4 doublezt = target->getZ();1148 G4int ab = bullet->getA(); 1149 G4int zb = bullet->getZ(); 1150 G4int at = target->getA(); 1151 G4int zt = target->getZ(); 1116 1152 1117 1153 G4double massb = bullet->getMass(); // For creating LorentzVectors below … … 1119 1155 if (ab < max_a_for_cascad) { 1120 1156 1121 G4double benb = 0.001 * G4NucleiProperties::GetBindingEnergy(G4lrint(ab), G4lrint(zb)) / ab;1122 G4double bent = 0.001 * G4NucleiProperties::GetBindingEnergy(G4lrint(at), G4lrint(zt)) / at;1157 G4double benb = 0.001 * bindingEnergy(ab,zb) / G4double(ab); 1158 G4double bent = 0.001 * bindingEnergy(at,zt) / G4double(at); 1123 1159 G4double ben = benb < bent ? bent : benb; 1124 1160 … … 1128 1164 while (casparticles.size() == 0 && itryg < itry_max) { 1129 1165 itryg++; 1130 1131 if(itryg > 0) particles.clear(); 1166 particles.clear(); 1132 1167 1133 1168 // nucleons coordinates and momenta in nuclei rest frame 1134 std::vector<G4ThreeVector> coordinates;1135 std::vector<G4LorentzVector> momentums;1169 coordinates.clear(); 1170 momentums.clear(); 1136 1171 1137 if (ab < 3 .0) { // deutron, simplest case1172 if (ab < 3) { // deuteron, simplest case 1138 1173 G4double r = 2.214 - 3.4208 * std::log(1.0 - 0.981 * inuclRndm()); 1139 G4double s = 2.0 * inuclRndm() - 1.0; 1140 G4double r1 = r * std::sqrt(1.0 - s * s); 1141 G4double phi = randomPHI(); 1142 1143 G4ThreeVector coord1(r1*std::cos(phi), r1*std::sin(phi), r*s); 1174 G4ThreeVector coord1 = generateWithRandomAngles(r).vect(); 1144 1175 coordinates.push_back(coord1); 1145 1146 coord1 *= -1.; 1147 coordinates.push_back(coord1); 1176 coordinates.push_back(-coord1); 1148 1177 1149 1178 G4double p = 0.0; … … 1157 1186 if (p * p / (p * p + 2079.36) / (p * p + 2079.36) > 1.2023e-4 * inuclRndm() && 1158 1187 p * r > 312.0) bad = false; 1159 } ;1188 } 1160 1189 1161 1190 if (itry == itry_max) … … 1176 1205 momentums.push_back(-mom); 1177 1206 } else { 1178 G4int ia = int(ab + 0.5);1179 1180 1207 G4ThreeVector coord1; 1181 1208 … … 1184 1211 G4int itry = 0; 1185 1212 1186 if (ab < 4.0) { // a == 31213 if (ab == 3) { 1187 1214 while (badco && itry < itry_max) { 1188 1215 if (itry > 0) coordinates.clear(); … … 1210 1237 } 1211 1238 break; 1212 } ;1213 } ;1239 } 1240 } 1214 1241 1215 1242 if (itry1 == itry_max) { // bad case … … 1217 1244 coordinates.push_back(coord1); 1218 1245 break; 1219 } ;1220 } ;1246 } 1247 } 1221 1248 1222 1249 coord1 = -coordinates[0] - coordinates[1]; … … 1241 1268 1242 1269 break; 1243 } ;1244 } ;1270 } 1271 } 1245 1272 1246 1273 if (large_dist) break; 1247 } ;1274 } 1248 1275 1249 1276 if(!large_dist) badco = false; 1250 1277 1251 } ;1278 } 1252 1279 1253 1280 } else { // a >= 4 … … 1264 1291 G4int i(0); 1265 1292 1266 for (i = 0; i < ia-1; i++) {1293 for (i = 0; i < ab-1; i++) { 1267 1294 G4int itry1 = 0; 1268 1295 G4double s, u; … … 1283 1310 1284 1311 break; 1285 } ;1286 } ;1312 } 1313 } 1287 1314 1288 1315 if (itry1 == itry_max) { // bad case … … 1290 1317 coordinates.push_back(coord1); 1291 1318 break; 1292 } ;1293 } ;1319 } 1320 } 1294 1321 1295 1322 coord1 *= 0.0; // Cheap way to reset 1296 for(G4int j = 0; j < ia-1; j++) coord1 -= coordinates[j];1323 for(G4int j = 0; j < ab -1; j++) coord1 -= coordinates[j]; 1297 1324 1298 1325 coordinates.push_back(coord1); … … 1304 1331 G4bool large_dist = false; 1305 1332 1306 for (i = 0; i < ia-1; i++) {1307 for (G4int j = i+1; j < ia; j++) {1333 for (i = 0; i < ab-1; i++) { 1334 for (G4int j = i+1; j < ab; j++) { 1308 1335 1309 1336 G4double r2 = (coordinates[i]-coordinates[j]).mag2(); … … 1317 1344 1318 1345 break; 1319 } ;1320 } ;1346 } 1347 } 1321 1348 1322 1349 if (large_dist) break; 1323 } ;1350 } 1324 1351 1325 1352 if (!large_dist) badco = false; 1326 } ;1327 } ;1353 } 1354 } 1328 1355 1329 1356 if(badco) { … … 1339 1366 G4LorentzVector mom; 1340 1367 //G4bool badp = True; 1341 G4int i(0); 1342 1343 for (i = 0; i < ia - 1; i++) { 1368 1369 for (G4int i = 0; i < ab - 1; i++) { 1344 1370 G4int itry = 0; 1345 1371 … … 1355 1381 1356 1382 break; 1357 } ;1358 } ;1383 } 1384 } 1359 1385 1360 1386 if(itry == itry_max) { … … 1365 1391 particles.clear(); 1366 1392 return; 1367 } ;1368 1369 } ;1393 } 1394 1395 } 1370 1396 // last momentum 1371 1397 1372 1398 mom *= 0.; // Cheap way to reset 1373 for(G4int j=0; j< ia-1; j++) mom -= momentums[j]; 1399 mom.setE(bullet->getEnergy()+target->getEnergy()); 1400 1401 for(G4int j=0; j< ab-1; j++) mom -= momentums[j]; 1374 1402 1375 1403 momentums.push_back(mom); 1376 } ;1404 } 1377 1405 } 1378 1406 … … 1385 1413 1386 1414 if(rp > rb) rb = rp; 1387 } ;1415 } 1388 1416 1389 1417 // nuclei i.p. as a whole … … 1396 1424 for (i = 0; i < G4int(coordinates.size()); i++) { 1397 1425 coordinates[i] += global_pos; 1398 } ;1426 } 1399 1427 1400 1428 // all nucleons at rest 1401 std::vector<G4InuclElementaryParticle> raw_particles; 1402 G4int ia = G4int(ab + 0.5); 1403 G4int iz = G4int(zb + 0.5); 1404 1405 for (G4int ipa = 0; ipa < ia; ipa++) { 1406 G4int knd = ipa < iz ? 1 : 2; 1429 raw_particles.clear(); 1430 1431 for (G4int ipa = 0; ipa < ab; ipa++) { 1432 G4int knd = ipa < zb ? 1 : 2; 1407 1433 raw_particles.push_back(G4InuclElementaryParticle(momentums[ipa], knd)); 1408 } ;1434 } 1409 1435 1410 1436 G4InuclElementaryParticle dummy(small_ekin, 1); 1411 G4LorentzConvertor toTheBulletRestFrame; 1412 toTheBulletRestFrame.setBullet(dummy.getMomentum(), dummy.getMass()); 1413 toTheBulletRestFrame.setTarget(bullet->getMomentum(),bullet->getMass()); 1437 G4LorentzConvertor toTheBulletRestFrame(&dummy, bullet); 1414 1438 toTheBulletRestFrame.toTheTargetRestFrame(); 1415 1439 … … 1418 1442 for (ipart = raw_particles.begin(); ipart != raw_particles.end(); ipart++) { 1419 1443 ipart->setMomentum(toTheBulletRestFrame.backToTheLab(ipart->getMomentum())); 1420 } ;1444 } 1421 1445 1422 1446 // fill cascad particles and outgoing particles … … 1437 1461 if(t1 > 0.0) { 1438 1462 if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1; 1439 } ;1463 } 1440 1464 1441 1465 if(tr < 0.0 && t2 > 0.0) { 1442 1466 1443 1467 if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2; 1444 } ;1468 } 1445 1469 1446 1470 } else { … … 1448 1472 1449 1473 if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2; 1450 } ;1474 } 1451 1475 1452 1476 if(tr < 0.0 && t1 > 0.0) { 1453 1477 if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1; 1454 } ;1455 } ;1456 1457 } ;1478 } 1479 } 1480 1481 } 1458 1482 1459 1483 if(tr >= 0.0) { // cascad particle 1460 coordinates[ip] += mom *tr / pmod;1484 coordinates[ip] += mom.vect()*tr / pmod; 1461 1485 casparticles.push_back(G4CascadParticle(raw_particles[ip], 1462 1486 coordinates[ip], … … 1465 1489 } else { 1466 1490 particles.push_back(raw_particles[ip]); 1467 } ;1468 } ;1469 } ;1491 } 1492 } 1493 } 1470 1494 1471 1495 if(casparticles.size() == 0) { … … 1474 1498 G4cout << " can not generate proper distribution for " << itry_max 1475 1499 << " steps " << G4endl; 1476 } ;1477 } ;1478 } ;1500 } 1501 } 1502 } 1479 1503 1480 1504 if(verboseLevel > 2){ … … 1509 1533 } else if (ke < 1.0) { 1510 1534 csec = 3.6735 * (1.0-ke)*(1.0-ke); 1511 } ;1535 } 1512 1536 1513 1537 if (csec < 0.0) csec = 0.0;
Note: See TracChangeset
for help on using the changeset viewer.