00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116 #include <numeric>
00117
00118 #include "G4NucleiModel.hh"
00119 #include "G4PhysicalConstants.hh"
00120 #include "G4SystemOfUnits.hh"
00121 #include "G4CascadeChannel.hh"
00122 #include "G4CascadeChannelTables.hh"
00123 #include "G4CascadeCheckBalance.hh"
00124 #include "G4CascadeInterpolator.hh"
00125 #include "G4CascadeParameters.hh"
00126 #include "G4CollisionOutput.hh"
00127 #include "G4ElementaryParticleCollider.hh"
00128 #include "G4HadTmpUtil.hh"
00129 #include "G4InuclNuclei.hh"
00130 #include "G4InuclParticleNames.hh"
00131 #include "G4InuclSpecialFunctions.hh"
00132 #include "G4LorentzConvertor.hh"
00133 #include "G4Neutron.hh"
00134 #include "G4ParticleLargerBeta.hh"
00135 #include "G4ParticleDefinition.hh"
00136 #include "G4Proton.hh"
00137
00138 using namespace G4InuclParticleNames;
00139 using namespace G4InuclSpecialFunctions;
00140
00141 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153 const G4double G4NucleiModel::crossSectionUnits = G4CascadeParameters::xsecScale();
00154 const G4double G4NucleiModel::radiusUnits = G4CascadeParameters::radiusScale();
00155 const G4double G4NucleiModel::skinDepth = 0.611207*radiusUnits;
00156
00157
00158
00159
00160 const G4double G4NucleiModel::radiusScale =
00161 (G4CascadeParameters::useTwoParam() ? 1.16 : 1.2) * radiusUnits;
00162 const G4double G4NucleiModel::radiusScale2 =
00163 (G4CascadeParameters::useTwoParam() ? -1.3456 : 0.) * radiusUnits;
00164
00165
00166
00167
00168 const G4double G4NucleiModel::radiusForSmall = G4CascadeParameters::radiusSmall();
00169
00170 const G4double G4NucleiModel::radScaleAlpha = G4CascadeParameters::radiusAlpha();
00171
00172
00173
00174
00175 const G4double G4NucleiModel::fermiMomentum = G4CascadeParameters::fermiScale();
00176
00177
00178 const G4double G4NucleiModel::R_nucleon = G4CascadeParameters::radiusTrailing();
00179
00180
00181 const G4double G4NucleiModel::alfa3[3] = { 0.7, 0.3, 0.01 };
00182 const G4double G4NucleiModel::alfa6[6] = { 0.9, 0.6, 0.4, 0.2, 0.1, 0.05 };
00183
00184
00185 const G4double G4NucleiModel::pion_vp = 0.007;
00186 const G4double G4NucleiModel::pion_vp_small = 0.007;
00187 const G4double G4NucleiModel::kaon_vp = 0.015;
00188 const G4double G4NucleiModel::hyperon_vp = 0.030;
00189
00190 const G4double G4NucleiModel::piTimes4thirds = pi*4./3.;
00191
00192
00193
00194 G4NucleiModel::G4NucleiModel()
00195 : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
00196 A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
00197 neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
00198 current_nucl2(0) {}
00199
00200 G4NucleiModel::G4NucleiModel(G4int a, G4int z)
00201 : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
00202 A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
00203 neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
00204 current_nucl2(0) {
00205 generateModel(a,z);
00206 }
00207
00208 G4NucleiModel::G4NucleiModel(G4InuclNuclei* nuclei)
00209 : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
00210 A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
00211 neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
00212 current_nucl2(0) {
00213 generateModel(nuclei);
00214 }
00215
00216 G4NucleiModel::~G4NucleiModel() {
00217 delete theNucleus;
00218 theNucleus = 0;
00219 }
00220
00221
00222
00223
00224 void G4NucleiModel::reset(G4int nHitNeutrons, G4int nHitProtons,
00225 const std::vector<G4ThreeVector>* hitPoints) {
00226 neutronNumberCurrent = neutronNumber - nHitNeutrons;
00227 protonNumberCurrent = protonNumber - nHitProtons;
00228
00229
00230 if (!hitPoints || !hitPoints->empty()) collisionPts.clear();
00231 else collisionPts = *hitPoints;
00232 }
00233
00234
00235
00236
00237 void G4NucleiModel::generateModel(G4InuclNuclei* nuclei) {
00238 generateModel(nuclei->getA(), nuclei->getZ());
00239 }
00240
00241 void G4NucleiModel::generateModel(G4int a, G4int z) {
00242 if (verboseLevel) {
00243 G4cout << " >>> G4NucleiModel::generateModel A " << a << " Z " << z
00244 << G4endl;
00245 }
00246
00247
00248 if (a == A && z == Z) {
00249 if (verboseLevel > 1) G4cout << " model already generated" << z << G4endl;
00250 reset();
00251 return;
00252 }
00253
00254 A = a;
00255 Z = z;
00256 delete theNucleus;
00257 theNucleus = new G4InuclNuclei(A,Z);
00258
00259 neutronNumber = A - Z;
00260 protonNumber = Z;
00261 reset();
00262
00263 if (verboseLevel > 3) {
00264 G4cout << " crossSectionUnits = " << crossSectionUnits << G4endl
00265 << " radiusUnits = " << radiusUnits << G4endl
00266 << " skinDepth = " << skinDepth << G4endl
00267 << " radiusScale = " << radiusScale << G4endl
00268 << " radiusScale2 = " << radiusScale2 << G4endl
00269 << " radiusForSmall = " << radiusForSmall << G4endl
00270 << " radScaleAlpha = " << radScaleAlpha << G4endl
00271 << " fermiMomentum = " << fermiMomentum << G4endl
00272 << " piTimes4thirds = " << piTimes4thirds << G4endl;
00273 }
00274
00275 G4double nuclearRadius;
00276 if (A>4) nuclearRadius = radiusScale*G4cbrt(A) + radiusScale2/G4cbrt(A);
00277 else nuclearRadius = radiusForSmall * (A==4 ? radScaleAlpha : 1.);
00278
00279
00280 number_of_zones = (A < 5) ? 1 : (A < 100) ? 3 : 6;
00281
00282
00283 binding_energies.clear();
00284 nucleon_densities.clear();
00285 zone_potentials.clear();
00286 fermi_momenta.clear();
00287 zone_radii.clear();
00288 zone_volumes.clear();
00289
00290 fillBindingEnergies();
00291 fillZoneRadii(nuclearRadius);
00292
00293 G4double tot_vol = fillZoneVolumes(nuclearRadius);
00294
00295 fillPotentials(proton, tot_vol);
00296 fillPotentials(neutron, tot_vol);
00297
00298
00299 const std::vector<G4double> vp(number_of_zones, (A>4)?pion_vp:pion_vp_small);
00300 const std::vector<G4double> kp(number_of_zones, kaon_vp);
00301 const std::vector<G4double> hp(number_of_zones, hyperon_vp);
00302
00303 zone_potentials.push_back(vp);
00304 zone_potentials.push_back(kp);
00305 zone_potentials.push_back(hp);
00306
00307 nuclei_radius = zone_radii.back();
00308 nuclei_volume = std::accumulate(zone_volumes.begin(),zone_volumes.end(),0.);
00309
00310 if (verboseLevel > 3) printModel();
00311 }
00312
00313
00314
00315
00316 void G4NucleiModel::fillBindingEnergies() {
00317 if (verboseLevel > 1)
00318 G4cout << " >>> G4NucleiModel::fillBindingEnergies" << G4endl;
00319
00320 G4double dm = bindingEnergy(A,Z);
00321
00322
00323
00324 binding_energies.push_back(std::fabs(bindingEnergy(A-1,Z-1)-dm)/GeV);
00325 binding_energies.push_back(std::fabs(bindingEnergy(A-1,Z)-dm)/GeV);
00326 }
00327
00328
00329
00330 void G4NucleiModel::fillZoneRadii(G4double nuclearRadius) {
00331 if (verboseLevel > 1)
00332 G4cout << " >>> G4NucleiModel::fillZoneRadii" << G4endl;
00333
00334 G4double skinRatio = nuclearRadius/skinDepth;
00335 G4double skinDecay = std::exp(-skinRatio);
00336
00337 if (A < 5) {
00338 zone_radii.push_back(nuclearRadius);
00339 ur[0] = 0.;
00340 ur[1] = 1.;
00341 } else if (A < 12) {
00342 G4double rSq = nuclearRadius * nuclearRadius;
00343 G4double gaussRadius = std::sqrt(rSq * (1.0 - 1.0/A) + 6.4);
00344
00345 ur[0] = 0.0;
00346 for (G4int i = 0; i < number_of_zones; i++) {
00347 G4double y = std::sqrt(-std::log(alfa3[i]));
00348 zone_radii.push_back(gaussRadius * y);
00349 ur[i+1] = y;
00350 }
00351 } else if (A < 100) {
00352 ur[0] = -skinRatio;
00353 for (G4int i = 0; i < number_of_zones; i++) {
00354 G4double y = std::log((1.0 + skinDecay)/alfa3[i] - 1.0);
00355 zone_radii.push_back(nuclearRadius + skinDepth * y);
00356 ur[i+1] = y;
00357 }
00358 } else {
00359 ur[0] = -skinRatio;
00360 for (G4int i = 0; i < number_of_zones; i++) {
00361 G4double y = std::log((1.0 + skinDecay)/alfa6[i] - 1.0);
00362 zone_radii.push_back(nuclearRadius + skinDepth * y);
00363 ur[i+1] = y;
00364 }
00365 }
00366 }
00367
00368
00369
00370 G4double G4NucleiModel::fillZoneVolumes(G4double nuclearRadius) {
00371 if (verboseLevel > 1)
00372 G4cout << " >>> G4NucleiModel::fillZoneVolumes" << G4endl;
00373
00374 G4double tot_vol = 0.;
00375
00376 if (A < 5) {
00377 v[0] = v1[0] = 1.;
00378 tot_vol = zone_radii[0]*zone_radii[0]*zone_radii[0];
00379 zone_volumes.push_back(tot_vol*piTimes4thirds);
00380
00381 return tot_vol;
00382 }
00383
00384 PotentialType usePotential = (A < 12) ? Gaussian : WoodsSaxon;
00385
00386 for (G4int i = 0; i < number_of_zones; i++) {
00387 if (usePotential == WoodsSaxon) {
00388 v[i] = zoneIntegralWoodsSaxon(ur[i], ur[i+1], nuclearRadius);
00389 } else {
00390 v[i] = zoneIntegralGaussian(ur[i], ur[i+1], nuclearRadius);
00391 }
00392
00393 tot_vol += v[i];
00394
00395 v1[i] = zone_radii[i]*zone_radii[i]*zone_radii[i];
00396 if (i > 0) v1[i] -= zone_radii[i-1]*zone_radii[i-1]*zone_radii[i-1];
00397
00398 zone_volumes.push_back(v1[i]*piTimes4thirds);
00399 }
00400
00401 return tot_vol;
00402 }
00403
00404
00405 void G4NucleiModel::fillPotentials(G4int type, G4double tot_vol) {
00406 if (verboseLevel > 1)
00407 G4cout << " >>> G4NucleiModel::fillZoneVolumes(" << type << ")" << G4endl;
00408
00409 if (type != proton && type != neutron) return;
00410
00411 const G4double mass = G4InuclElementaryParticle::getParticleMass(type);
00412
00413
00414 const G4double dm = binding_energies[type-1];
00415
00416 rod.clear(); rod.reserve(number_of_zones);
00417 pf.clear(); pf.reserve(number_of_zones);
00418 vz.clear(); vz.reserve(number_of_zones);
00419
00420 G4int nNucleons = (type==proton) ? protonNumber : neutronNumber;
00421 G4double dd0 = nNucleons / tot_vol / piTimes4thirds;
00422
00423 for (G4int i = 0; i < number_of_zones; i++) {
00424 G4double rd = dd0 * v[i] / v1[i];
00425 rod.push_back(rd);
00426 G4double pff = fermiMomentum * G4cbrt(rd);
00427 pf.push_back(pff);
00428 vz.push_back(0.5 * pff * pff / mass + dm);
00429 }
00430
00431 nucleon_densities.push_back(rod);
00432 fermi_momenta.push_back(pf);
00433 zone_potentials.push_back(vz);
00434 }
00435
00436
00437 G4double G4NucleiModel::zoneIntegralWoodsSaxon(G4double r1, G4double r2,
00438 G4double nuclearRadius) const {
00439 if (verboseLevel > 1) {
00440 G4cout << " >>> G4NucleiModel::zoneIntegralWoodsSaxon" << G4endl;
00441 }
00442
00443 const G4double epsilon = 1.0e-3;
00444 const G4int itry_max = 1000;
00445
00446 G4double skinRatio = nuclearRadius / skinDepth;
00447
00448 G4double d2 = 2.0 * skinRatio;
00449 G4double dr = r2 - r1;
00450 G4double fr1 = r1 * (r1 + d2) / (1.0 + std::exp(r1));
00451 G4double fr2 = r2 * (r2 + d2) / (1.0 + std::exp(r2));
00452 G4double fi = (fr1 + fr2) / 2.;
00453 G4double fun1 = fi * dr;
00454 G4double fun;
00455 G4double jc = 1;
00456 G4double dr1 = dr;
00457 G4int itry = 0;
00458
00459 while (itry < itry_max) {
00460 dr /= 2.;
00461 itry++;
00462
00463 G4double r = r1 - dr;
00464 fi = 0.0;
00465 G4int jc1 = G4int(std::pow(2.0, jc - 1) + 0.1);
00466
00467 for (G4int i = 0; i < jc1; i++) {
00468 r += dr1;
00469 fi += r * (r + d2) / (1.0 + std::exp(r));
00470 }
00471
00472 fun = 0.5 * fun1 + fi * dr;
00473
00474 if (std::fabs((fun - fun1) / fun) <= epsilon) break;
00475
00476 jc++;
00477 dr1 = dr;
00478 fun1 = fun;
00479 }
00480
00481 if (verboseLevel > 2 && itry == itry_max)
00482 G4cout << " zoneIntegralWoodsSaxon-> n iter " << itry_max << G4endl;
00483
00484 G4double skinDepth3 = skinDepth*skinDepth*skinDepth;
00485
00486 return skinDepth3 * (fun + skinRatio*skinRatio*std::log((1.0 + std::exp(-r1)) / (1.0 + std::exp(-r2))));
00487 }
00488
00489
00490
00491 G4double G4NucleiModel::zoneIntegralGaussian(G4double r1, G4double r2,
00492 G4double nucRad) const {
00493 if (verboseLevel > 1) {
00494 G4cout << " >>> G4NucleiModel::zoneIntegralGaussian" << G4endl;
00495 }
00496
00497 G4double gaussRadius = std::sqrt(nucRad*nucRad * (1.0 - 1.0/A) + 6.4);
00498
00499 const G4double epsilon = 1.0e-3;
00500 const G4int itry_max = 1000;
00501
00502 G4double dr = r2 - r1;
00503 G4double fr1 = r1 * r1 * std::exp(-r1 * r1);
00504 G4double fr2 = r2 * r2 * std::exp(-r2 * r2);
00505 G4double fi = (fr1 + fr2) / 2.;
00506 G4double fun1 = fi * dr;
00507 G4double fun;
00508 G4double jc = 1;
00509 G4double dr1 = dr;
00510 G4int itry = 0;
00511
00512 while (itry < itry_max) {
00513 dr /= 2.;
00514 itry++;
00515 G4double r = r1 - dr;
00516 fi = 0.0;
00517 G4int jc1 = int(std::pow(2.0, jc - 1) + 0.1);
00518
00519 for (G4int i = 0; i < jc1; i++) {
00520 r += dr1;
00521 fi += r * r * std::exp(-r * r);
00522 }
00523
00524 fun = 0.5 * fun1 + fi * dr;
00525
00526 if (std::fabs((fun - fun1) / fun) <= epsilon) break;
00527
00528 jc++;
00529 dr1 = dr;
00530 fun1 = fun;
00531 }
00532
00533 if (verboseLevel > 2 && itry == itry_max)
00534 G4cerr << " zoneIntegralGaussian-> n iter " << itry_max << G4endl;
00535
00536 return gaussRadius*gaussRadius*gaussRadius * fun;
00537 }
00538
00539
00540 void G4NucleiModel::printModel() const {
00541 if (verboseLevel > 1) {
00542 G4cout << " >>> G4NucleiModel::printModel" << G4endl;
00543 }
00544
00545 G4cout << " nuclei model for A " << A << " Z " << Z << G4endl
00546 << " proton binding energy " << binding_energies[0]
00547 << " neutron binding energy " << binding_energies[1] << G4endl
00548 << " Nuclei radius " << nuclei_radius << " volume " << nuclei_volume
00549 << " number of zones " << number_of_zones << G4endl;
00550
00551 for (G4int i = 0; i < number_of_zones; i++)
00552 G4cout << " zone " << i+1 << " radius " << zone_radii[i]
00553 << " volume " << zone_volumes[i] << G4endl
00554 << " protons: density " << getDensity(1,i) << " PF " <<
00555 getFermiMomentum(1,i) << " VP " << getPotential(1,i) << G4endl
00556 << " neutrons: density " << getDensity(2,i) << " PF " <<
00557 getFermiMomentum(2,i) << " VP " << getPotential(2,i) << G4endl
00558 << " pions: VP " << getPotential(3,i) << G4endl;
00559 }
00560
00561
00562 G4double G4NucleiModel::getFermiKinetic(G4int ip, G4int izone) const {
00563 G4double ekin = 0.0;
00564
00565 if (ip < 3 && izone < number_of_zones) {
00566 G4double pfermi = fermi_momenta[ip - 1][izone];
00567 G4double mass = G4InuclElementaryParticle::getParticleMass(ip);
00568
00569 ekin = std::sqrt(pfermi*pfermi + mass*mass) - mass;
00570 }
00571 return ekin;
00572 }
00573
00574
00575 G4LorentzVector
00576 G4NucleiModel::generateNucleonMomentum(G4int type, G4int zone) const {
00577 G4double pmod = getFermiMomentum(type, zone) * G4cbrt(inuclRndm());
00578 G4double mass = G4InuclElementaryParticle::getParticleMass(type);
00579
00580 return generateWithRandomAngles(pmod, mass);
00581 }
00582
00583
00584 G4InuclElementaryParticle
00585 G4NucleiModel::generateNucleon(G4int type, G4int zone) const {
00586 if (verboseLevel > 1) {
00587 G4cout << " >>> G4NucleiModel::generateNucleon" << G4endl;
00588 }
00589
00590 G4LorentzVector mom = generateNucleonMomentum(type, zone);
00591 return G4InuclElementaryParticle(mom, type);
00592 }
00593
00594
00595 G4InuclElementaryParticle
00596 G4NucleiModel::generateQuasiDeutron(G4int type1, G4int type2,
00597 G4int zone) const {
00598 if (verboseLevel > 1) {
00599 G4cout << " >>> G4NucleiModel::generateQuasiDeutron" << G4endl;
00600 }
00601
00602
00603
00604
00605
00606
00607
00608 G4LorentzVector mom1 = generateNucleonMomentum(type1, zone);
00609 G4LorentzVector mom2 = generateNucleonMomentum(type2, zone);
00610 G4LorentzVector dmom = mom1+mom2;
00611
00612 G4int dtype = 0;
00613 if (type1*type2 == pro*pro) dtype = 111;
00614 else if (type1*type2 == pro*neu) dtype = 112;
00615 else if (type1*type2 == neu*neu) dtype = 122;
00616
00617 return G4InuclElementaryParticle(dmom, dtype);
00618 }
00619
00620
00621 void
00622 G4NucleiModel::generateInteractionPartners(G4CascadParticle& cparticle) {
00623 if (verboseLevel > 1) {
00624 G4cout << " >>> G4NucleiModel::generateInteractionPartners" << G4endl;
00625 }
00626
00627 const G4double small = 1.0e-10;
00628
00629 thePartners.clear();
00630
00631 G4int ptype = cparticle.getParticle().type();
00632 G4int zone = cparticle.getCurrentZone();
00633
00634 G4double r_in;
00635 G4double r_out;
00636
00637 if (zone == number_of_zones) {
00638 r_in = nuclei_radius;
00639 r_out = 0.0;
00640 } else if (zone == 0) {
00641 r_in = 0.0;
00642 r_out = zone_radii[0];
00643 } else {
00644 r_in = zone_radii[zone - 1];
00645 r_out = zone_radii[zone];
00646 }
00647
00648 G4double path = cparticle.getPathToTheNextZone(r_in, r_out);
00649
00650 if (verboseLevel > 2) {
00651 if (isProjectile(cparticle)) G4cout << " incident particle" << G4endl;
00652 G4cout << " r_in " << r_in << " r_out " << r_out << " path " << path
00653 << G4endl;
00654 }
00655
00656 if (path < -small) {
00657 if (verboseLevel)
00658 G4cerr << " generateInteractionPartners-> negative path length" << G4endl;
00659 return;
00660 }
00661
00662 if (std::fabs(path) < small) {
00663 if (verboseLevel > 3)
00664 G4cout << " generateInteractionPartners-> zero path" << G4endl;
00665
00666 thePartners.push_back(partner());
00667 return;
00668 }
00669
00670
00671 dummy_convertor.setBullet(cparticle.getParticle());
00672
00673 G4double invmfp = 0.;
00674 G4double spath = 0.;
00675 for (G4int ip = 1; ip < 3; ip++) {
00676
00677 if (ip==proton && protonNumberCurrent < 1) continue;
00678 if (ip==neutron && neutronNumberCurrent < 1) continue;
00679
00680
00681 G4InuclElementaryParticle particle = generateNucleon(ip, zone);
00682 invmfp = inverseMeanFreePath(cparticle, particle);
00683 spath = generateInteractionLength(cparticle, path, invmfp);
00684
00685 if (spath < path) {
00686 if (verboseLevel > 3) {
00687 G4cout << " adding partner[" << thePartners.size() << "]: "
00688 << particle << G4endl;
00689 }
00690 thePartners.push_back(partner(particle, spath));
00691 }
00692 }
00693
00694 if (verboseLevel > 2) {
00695 G4cout << " after nucleons " << thePartners.size() << " path " << path << G4endl;
00696 }
00697
00698
00699 if (cparticle.getParticle().pion() || cparticle.getParticle().isPhoton()) {
00700 if (verboseLevel > 2) {
00701 G4cout << " trying quasi-deuterons with bullet: "
00702 << cparticle.getParticle() << G4endl;
00703 }
00704
00705
00706 qdeutrons.clear();
00707 acsecs.clear();
00708
00709 G4double tot_invmfp = 0.0;
00710
00711
00712 if (protonNumberCurrent >= 2 && ptype != pip) {
00713 G4InuclElementaryParticle ppd = generateQuasiDeutron(pro, pro, zone);
00714 if (verboseLevel > 2)
00715 G4cout << " ptype=" << ptype << " using pp target\n" << ppd << G4endl;
00716
00717 invmfp = inverseMeanFreePath(cparticle, ppd);
00718 tot_invmfp += invmfp;
00719 acsecs.push_back(invmfp);
00720 qdeutrons.push_back(ppd);
00721 }
00722
00723
00724 if (protonNumberCurrent >= 1 && neutronNumberCurrent >= 1) {
00725 G4InuclElementaryParticle npd = generateQuasiDeutron(pro, neu, zone);
00726 if (verboseLevel > 2)
00727 G4cout << " ptype=" << ptype << " using np target\n" << npd << G4endl;
00728
00729 invmfp = inverseMeanFreePath(cparticle, npd);
00730 tot_invmfp += invmfp;
00731 acsecs.push_back(invmfp);
00732 qdeutrons.push_back(npd);
00733 }
00734
00735
00736 if (neutronNumberCurrent >= 2 && ptype != pim) {
00737 G4InuclElementaryParticle nnd = generateQuasiDeutron(neu, neu, zone);
00738 if (verboseLevel > 2)
00739 G4cout << " ptype=" << ptype << " using nn target\n" << nnd << G4endl;
00740
00741 invmfp = inverseMeanFreePath(cparticle, nnd);
00742 tot_invmfp += invmfp;
00743 acsecs.push_back(invmfp);
00744 qdeutrons.push_back(nnd);
00745 }
00746
00747
00748 if (verboseLevel > 2) {
00749 for (size_t i=0; i<qdeutrons.size(); i++) {
00750 G4cout << " acsecs[" << qdeutrons[i].getDefinition()->GetParticleName()
00751 << "] " << acsecs[i];
00752 }
00753 G4cout << G4endl;
00754 }
00755
00756 if (tot_invmfp > small) {
00757 G4double apath = generateInteractionLength(cparticle, path, tot_invmfp);
00758
00759 if (apath < path) {
00760 G4double sl = inuclRndm() * tot_invmfp;
00761 G4double as = 0.0;
00762
00763 for (size_t i = 0; i < qdeutrons.size(); i++) {
00764 as += acsecs[i];
00765 if (sl < as) {
00766 if (verboseLevel > 2) G4cout << " deut type " << i << G4endl;
00767 thePartners.push_back(partner(qdeutrons[i], apath));
00768 break;
00769 }
00770 }
00771 }
00772 }
00773 }
00774
00775 if (verboseLevel > 2) {
00776 G4cout << " after deuterons " << thePartners.size() << " partners"
00777 << G4endl;
00778 }
00779
00780 if (thePartners.size() > 1) {
00781 std::sort(thePartners.begin(), thePartners.end(), sortPartners);
00782 }
00783
00784 G4InuclElementaryParticle particle;
00785 thePartners.push_back(partner(particle, path));
00786 }
00787
00788
00789 void G4NucleiModel::
00790 generateParticleFate(G4CascadParticle& cparticle,
00791 G4ElementaryParticleCollider* theEPCollider,
00792 std::vector<G4CascadParticle>& outgoing_cparticles) {
00793 if (verboseLevel > 1)
00794 G4cout << " >>> G4NucleiModel::generateParticleFate" << G4endl;
00795
00796 if (verboseLevel > 2) G4cout << " cparticle: " << cparticle << G4endl;
00797
00798
00799 #ifdef G4CASCADE_CHECK_ECONS
00800 G4CascadeCheckBalance balance(0.005, 0.01, "G4NucleiModel");
00801 balance.setVerboseLevel(verboseLevel);
00802 #endif
00803
00804 outgoing_cparticles.clear();
00805 generateInteractionPartners(cparticle);
00806
00807 if (thePartners.empty()) {
00808 if (verboseLevel)
00809 G4cerr << " generateParticleFate-> got empty interaction-partners list "
00810 << G4endl;
00811 return;
00812 }
00813
00814 G4int npart = thePartners.size();
00815
00816 if (npart == 1) {
00817 cparticle.propagateAlongThePath(thePartners[0].second);
00818 cparticle.incrementCurrentPath(thePartners[0].second);
00819 boundaryTransition(cparticle);
00820 outgoing_cparticles.push_back(cparticle);
00821
00822 if (verboseLevel > 2) G4cout << " next zone \n" << cparticle << G4endl;
00823 } else {
00824 if (verboseLevel > 1)
00825 G4cout << " processing " << npart-1 << " possible interactions" << G4endl;
00826
00827 G4ThreeVector old_position = cparticle.getPosition();
00828 G4InuclElementaryParticle& bullet = cparticle.getParticle();
00829 G4bool no_interaction = true;
00830 G4int zone = cparticle.getCurrentZone();
00831
00832 for (G4int i=0; i<npart-1; i++) {
00833 if (i > 0) cparticle.updatePosition(old_position);
00834
00835 G4InuclElementaryParticle& target = thePartners[i].first;
00836
00837 if (verboseLevel > 3) {
00838 if (target.quasi_deutron()) G4cout << " try absorption: ";
00839 G4cout << " target " << target.type() << " bullet " << bullet.type()
00840 << G4endl;
00841 }
00842
00843 EPCoutput.reset();
00844 theEPCollider->collide(&bullet, &target, EPCoutput);
00845
00846
00847 if (EPCoutput.numberOfOutgoingParticles() == 0) break;
00848
00849 if (verboseLevel > 2) {
00850 EPCoutput.printCollisionOutput();
00851 #ifdef G4CASCADE_CHECK_ECONS
00852 balance.collide(&bullet, &target, EPCoutput);
00853 balance.okay();
00854 #endif
00855 }
00856
00857
00858 std::vector<G4InuclElementaryParticle>& outgoing_particles =
00859 EPCoutput.getOutgoingParticles();
00860
00861 if (!passFermi(outgoing_particles, zone)) continue;
00862
00863
00864 cparticle.propagateAlongThePath(thePartners[i].second);
00865 G4ThreeVector new_position = cparticle.getPosition();
00866
00867 if (!passTrailing(new_position)) continue;
00868 collisionPts.push_back(new_position);
00869
00870
00871 std::sort(outgoing_particles.begin(), outgoing_particles.end(),
00872 G4ParticleLargerBeta() );
00873
00874 if (verboseLevel > 2)
00875 G4cout << " adding " << outgoing_particles.size()
00876 << " output particles" << G4endl;
00877
00878
00879 for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++) {
00880 outgoing_cparticles.push_back(G4CascadParticle(outgoing_particles[ip],
00881 new_position, zone,
00882 0.0, 0));
00883 }
00884
00885 no_interaction = false;
00886 current_nucl1 = 0;
00887 current_nucl2 = 0;
00888
00889 #ifdef G4CASCADE_DEBUG_CHARGE
00890 {
00891 G4double out_charge = 0.0;
00892
00893 for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++)
00894 out_charge += outgoing_particles[ip].getCharge();
00895
00896 G4cout << " multiplicity " << outgoing_particles.size()
00897 << " bul type " << bullet.type()
00898 << " targ type " << target.type()
00899 << "\n initial charge "
00900 << bullet.getCharge() + target.getCharge()
00901 << " out charge " << out_charge << G4endl;
00902 }
00903 #endif
00904
00905 if (verboseLevel > 2)
00906 G4cout << " partner type " << target.type() << G4endl;
00907
00908 if (target.nucleon()) {
00909 current_nucl1 = target.type();
00910 } else {
00911 if (verboseLevel > 2) G4cout << " good absorption " << G4endl;
00912 current_nucl1 = (target.type() - 100) / 10;
00913 current_nucl2 = target.type() - 100 - 10 * current_nucl1;
00914 }
00915
00916 if (current_nucl1 == 1) {
00917 if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl;
00918 protonNumberCurrent--;
00919 } else {
00920 if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl;
00921 neutronNumberCurrent--;
00922 }
00923
00924 if (current_nucl2 == 1) {
00925 if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl;
00926 protonNumberCurrent--;
00927 } else if (current_nucl2 == 2) {
00928 if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl;
00929 neutronNumberCurrent--;
00930 }
00931
00932 break;
00933 }
00934
00935 if (no_interaction) {
00936 if (verboseLevel > 1) G4cout << " no interaction " << G4endl;
00937
00938
00939 static G4InuclElementaryParticle prescatCP;
00940 prescatCP = cparticle.getParticle();
00941
00942
00943 cparticle.updatePosition(old_position);
00944 cparticle.propagateAlongThePath(thePartners[npart-1].second);
00945 cparticle.incrementCurrentPath(thePartners[npart-1].second);
00946 boundaryTransition(cparticle);
00947 outgoing_cparticles.push_back(cparticle);
00948
00949
00950 #ifdef G4CASCADE_CHECK_ECONS
00951 if (verboseLevel > 2) {
00952 balance.collide(&prescatCP, 0, outgoing_cparticles);
00953 balance.okay();
00954 }
00955 #endif
00956 }
00957 }
00958
00959 return;
00960 }
00961
00962 G4bool G4NucleiModel::passFermi(const std::vector<G4InuclElementaryParticle>& particles,
00963 G4int zone) {
00964 if (verboseLevel > 1) {
00965 G4cout << " >>> G4NucleiModel::passFermi" << G4endl;
00966 }
00967
00968
00969 for (G4int i = 0; i < G4int(particles.size()); i++) {
00970 if (!particles[i].nucleon()) continue;
00971
00972 G4int type = particles[i].type();
00973 G4double mom = particles[i].getMomModule();
00974 G4double pfermi = fermi_momenta[type-1][zone];
00975
00976 if (verboseLevel > 2)
00977 G4cout << " type " << type << " p " << mom << " pf " << pfermi << G4endl;
00978
00979 if (mom < pfermi) {
00980 if (verboseLevel > 2) G4cout << " rejected by Fermi" << G4endl;
00981 return false;
00982 }
00983 }
00984 return true;
00985 }
00986
00987
00988
00989
00990
00991 G4bool G4NucleiModel::passTrailing(const G4ThreeVector& hit_position) {
00992 if (verboseLevel > 1)
00993 G4cout << " >>> G4NucleiModel::passTrailing " << hit_position << G4endl;
00994
00995 G4double dist;
00996 for (G4int i = 0; i < G4int(collisionPts.size() ); i++) {
00997 dist = (collisionPts[i] - hit_position).mag();
00998 if (verboseLevel > 2) G4cout << " dist " << dist << G4endl;
00999 if (dist < R_nucleon) {
01000 if (verboseLevel > 2) G4cout << " rejected by Trailing" << G4endl;
01001 return false;
01002 }
01003 }
01004 return true;
01005 }
01006
01007
01008 void G4NucleiModel::boundaryTransition(G4CascadParticle& cparticle) {
01009 if (verboseLevel > 1) {
01010 G4cout << " >>> G4NucleiModel::boundaryTransition" << G4endl;
01011 }
01012
01013 G4int zone = cparticle.getCurrentZone();
01014
01015 if (cparticle.movingInsideNuclei() && zone == 0) {
01016 G4cerr << " boundaryTransition-> in zone 0 " << G4endl;
01017 return;
01018 }
01019
01020 G4LorentzVector mom = cparticle.getMomentum();
01021 G4ThreeVector pos = cparticle.getPosition();
01022
01023 G4int type = cparticle.getParticle().type();
01024
01025 G4double r = pos.mag();
01026 G4double pr = pos.dot(mom.vect()) / r;
01027
01028 G4int next_zone = cparticle.movingInsideNuclei() ? zone - 1 : zone + 1;
01029
01030 G4double dv = getPotential(type,zone) - getPotential(type, next_zone);
01031 if (verboseLevel > 3) {
01032 G4cout << "Potentials for type " << type << " = "
01033 << getPotential(type,zone) << " , "
01034 << getPotential(type,next_zone) << G4endl;
01035 }
01036
01037 G4double qv = dv * dv - 2.0 * dv * mom.e() + pr * pr;
01038
01039 G4double p1r;
01040
01041 if (verboseLevel > 3) {
01042 G4cout << " type " << type << " zone " << zone << " next " << next_zone
01043 << " qv " << qv << " dv " << dv << G4endl;
01044 }
01045
01046 if (qv <= 0.0) {
01047 if (verboseLevel > 3) G4cout << " reflects off boundary" << G4endl;
01048 p1r = -pr;
01049 cparticle.incrementReflectionCounter();
01050 } else {
01051 if (verboseLevel > 3) G4cout << " passes thru boundary" << G4endl;
01052 p1r = std::sqrt(qv);
01053 if(pr < 0.0) p1r = -p1r;
01054 cparticle.updateZone(next_zone);
01055 cparticle.resetReflection();
01056 }
01057
01058 G4double prr = (p1r - pr) / r;
01059
01060 if (verboseLevel > 3) {
01061 G4cout << " prr " << prr << " delta px " << prr*pos.x() << " py "
01062 << prr*pos.y() << " pz " << prr*pos.z() << " mag "
01063 << std::fabs(prr*r) << G4endl;
01064 }
01065
01066 mom.setVect(mom.vect() + pos*prr);
01067 cparticle.updateParticleMomentum(mom);
01068 }
01069
01070
01071 G4bool G4NucleiModel::isProjectile(const G4CascadParticle& cparticle) const {
01072 if (verboseLevel > 2) {
01073 G4cout << " isProjectile(cparticle): zone "
01074 << cparticle.getCurrentZone() << " of " << number_of_zones
01075 << " path " << cparticle.getCurrentPath()
01076 << " movingInside " << cparticle.movingInsideNuclei()
01077 << " nReflections " << cparticle.getNumberOfReflections()
01078 << G4endl;
01079 }
01080
01081 return (cparticle.getCurrentZone() == number_of_zones-1 &&
01082 cparticle.getCurrentPath() == 1000. &&
01083 cparticle.getNumberOfReflections() <= 0 &&
01084 (cparticle.movingInsideNuclei()||number_of_zones==1)
01085 );
01086 }
01087
01088 G4bool G4NucleiModel::worthToPropagate(const G4CascadParticle& cparticle) const {
01089 if (verboseLevel > 1) {
01090 G4cout << " >>> G4NucleiModel::worthToPropagate" << G4endl;
01091 }
01092
01093 const G4double ekin_scale = 2.0;
01094
01095 G4bool worth = true;
01096
01097 if (cparticle.reflectedNow()) {
01098 G4int zone = cparticle.getCurrentZone();
01099 G4int ip = cparticle.getParticle().type();
01100
01101
01102 G4double ekin_cut = (cparticle.getParticle().nucleon()) ?
01103 getFermiKinetic(ip, zone) : 0.;
01104
01105 worth = cparticle.getParticle().getKineticEnergy()/ekin_scale > ekin_cut;
01106
01107 if (verboseLevel > 3) {
01108 G4cout << " type=" << ip
01109 << " ekin=" << cparticle.getParticle().getKineticEnergy()
01110 << " potential=" << ekin_cut
01111 << " : worth? " << worth << G4endl;
01112 }
01113 }
01114
01115 return worth;
01116 }
01117
01118
01119 G4double G4NucleiModel::getRatio(G4int ip) const {
01120 if (verboseLevel > 3) {
01121 G4cout << " >>> G4NucleiModel::getRatio " << ip << G4endl;
01122 }
01123
01124 switch (ip) {
01125 case proton: return G4double(protonNumberCurrent)/G4double(protonNumber);
01126 case neutron: return G4double(neutronNumberCurrent)/G4double(neutronNumber);
01127 case diproton: return getRatio(proton)*getRatio(proton);
01128 case unboundPN: return getRatio(proton)*getRatio(neutron);
01129 case dineutron: return getRatio(neutron)*getRatio(neutron);
01130 default: return 0.;
01131 }
01132
01133 return 0.;
01134 }
01135
01136 G4double G4NucleiModel::getCurrentDensity(G4int ip, G4int izone) const {
01137 const G4double pn_spec = 1.0;
01138
01139
01140 G4double dens = 0.;
01141
01142 if (ip < 100) dens = getDensity(ip,izone);
01143 else {
01144 switch (ip) {
01145 case diproton:
01146 dens = getDensity(proton,izone) * getDensity(proton,izone);
01147 break;
01148 case unboundPN:
01149 dens = getDensity(proton,izone) * getDensity(neutron,izone) * pn_spec;
01150 break;
01151 case dineutron:
01152 dens = getDensity(neutron,izone) * getDensity(neutron,izone);
01153 break;
01154 default: dens = 0.;
01155 }
01156 dens *= getVolume(izone);
01157 }
01158
01159 return getRatio(ip) * dens;
01160 }
01161
01162
01163 G4CascadParticle
01164 G4NucleiModel::initializeCascad(G4InuclElementaryParticle* particle) {
01165 if (verboseLevel > 1) {
01166 G4cout << " >>> G4NucleiModel::initializeCascad(particle)" << G4endl;
01167 }
01168
01169 const G4double large = 1000.0;
01170
01171
01172
01173
01174 G4double costh = std::sqrt(1.0 - inuclRndm());
01175 G4ThreeVector pos = generateWithFixedTheta(-costh, nuclei_radius);
01176
01177 G4CascadParticle cpart(*particle, pos, number_of_zones, large, 0);
01178
01179 if (verboseLevel > 2) G4cout << cpart << G4endl;
01180
01181 return cpart;
01182 }
01183
01184 void G4NucleiModel::initializeCascad(G4InuclNuclei* bullet,
01185 G4InuclNuclei* target,
01186 modelLists& output) {
01187 if (verboseLevel) {
01188 G4cout << " >>> G4NucleiModel::initializeCascad(bullet,target,output)"
01189 << G4endl;
01190 }
01191
01192 const G4double large = 1000.0;
01193 const G4double max_a_for_cascad = 5.0;
01194 const G4double ekin_cut = 2.0;
01195 const G4double small_ekin = 1.0e-6;
01196 const G4double r_large2for3 = 62.0;
01197 const G4double r0forAeq3 = 3.92;
01198 const G4double s3max = 6.5;
01199 const G4double r_large2for4 = 69.14;
01200 const G4double r0forAeq4 = 4.16;
01201 const G4double s4max = 7.0;
01202 const G4int itry_max = 100;
01203
01204
01205 std::vector<G4CascadParticle>& casparticles = output.first;
01206 std::vector<G4InuclElementaryParticle>& particles = output.second;
01207
01208 casparticles.clear();
01209 particles.clear();
01210
01211
01212 G4int ab = bullet->getA();
01213 G4int zb = bullet->getZ();
01214 G4int at = target->getA();
01215 G4int zt = target->getZ();
01216
01217 G4double massb = bullet->getMass();
01218
01219 if (ab < max_a_for_cascad) {
01220
01221 G4double benb = bindingEnergy(ab,zb)/GeV / G4double(ab);
01222 G4double bent = bindingEnergy(at,zt)/GeV / G4double(at);
01223 G4double ben = benb < bent ? bent : benb;
01224
01225 if (bullet->getKineticEnergy()/ab > ekin_cut*ben) {
01226 G4int itryg = 0;
01227
01228 while (casparticles.size() == 0 && itryg < itry_max) {
01229 itryg++;
01230 particles.clear();
01231
01232
01233 coordinates.clear();
01234 momentums.clear();
01235
01236 if (ab < 3) {
01237 G4double r = 2.214 - 3.4208 * std::log(1.0 - 0.981 * inuclRndm());
01238 G4ThreeVector coord1 = generateWithRandomAngles(r).vect();
01239 coordinates.push_back(coord1);
01240 coordinates.push_back(-coord1);
01241
01242 G4double p = 0.0;
01243 G4bool bad = true;
01244 G4int itry = 0;
01245
01246 while (bad && itry < itry_max) {
01247 itry++;
01248 p = 456.0 * inuclRndm();
01249
01250 if (p * p / (p * p + 2079.36) / (p * p + 2079.36) > 1.2023e-4 * inuclRndm() &&
01251 p * r > 312.0) bad = false;
01252 }
01253
01254 if (itry == itry_max)
01255 if (verboseLevel > 2){
01256 G4cout << " deutron bullet generation-> itry = " << itry_max << G4endl;
01257 }
01258
01259 p = 0.0005 * p;
01260
01261 if (verboseLevel > 2){
01262 G4cout << " p nuc " << p << G4endl;
01263 }
01264
01265 G4LorentzVector mom = generateWithRandomAngles(p, massb);
01266
01267 momentums.push_back(mom);
01268 mom.setVect(-mom.vect());
01269 momentums.push_back(-mom);
01270 } else {
01271 G4ThreeVector coord1;
01272
01273 G4bool badco = true;
01274
01275 G4int itry = 0;
01276
01277 if (ab == 3) {
01278 while (badco && itry < itry_max) {
01279 if (itry > 0) coordinates.clear();
01280 itry++;
01281 G4int i(0);
01282
01283 for (i = 0; i < 2; i++) {
01284 G4int itry1 = 0;
01285 G4double ss, u, rho;
01286 G4double fmax = std::exp(-0.5) / std::sqrt(0.5);
01287
01288 while (itry1 < itry_max) {
01289 itry1++;
01290 ss = -std::log(inuclRndm());
01291 u = fmax * inuclRndm();
01292 rho = std::sqrt(ss) * std::exp(-ss);
01293
01294 if (rho > u && ss < s3max) {
01295 ss = r0forAeq3 * std::sqrt(ss);
01296 coord1 = generateWithRandomAngles(ss).vect();
01297 coordinates.push_back(coord1);
01298
01299 if (verboseLevel > 2){
01300 G4cout << " i " << i << " r " << coord1.mag() << G4endl;
01301 }
01302 break;
01303 }
01304 }
01305
01306 if (itry1 == itry_max) {
01307 coord1.set(10000.,10000.,10000.);
01308 coordinates.push_back(coord1);
01309 break;
01310 }
01311 }
01312
01313 coord1 = -coordinates[0] - coordinates[1];
01314 if (verboseLevel > 2) {
01315 G4cout << " 3 r " << coord1.mag() << G4endl;
01316 }
01317
01318 coordinates.push_back(coord1);
01319
01320 G4bool large_dist = false;
01321
01322 for (i = 0; i < 2; i++) {
01323 for (G4int j = i+1; j < 3; j++) {
01324 G4double r2 = (coordinates[i]-coordinates[j]).mag2();
01325
01326 if (verboseLevel > 2) {
01327 G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
01328 }
01329
01330 if (r2 > r_large2for3) {
01331 large_dist = true;
01332
01333 break;
01334 }
01335 }
01336
01337 if (large_dist) break;
01338 }
01339
01340 if(!large_dist) badco = false;
01341
01342 }
01343
01344 } else {
01345 G4double b = 3./(ab - 2.0);
01346 G4double b1 = 1.0 - b / 2.0;
01347 G4double u = b1 + std::sqrt(b1 * b1 + b);
01348 G4double fmax = (1.0 + u/b) * u * std::exp(-u);
01349
01350 while (badco && itry < itry_max) {
01351
01352 if (itry > 0) coordinates.clear();
01353 itry++;
01354 G4int i(0);
01355
01356 for (i = 0; i < ab-1; i++) {
01357 G4int itry1 = 0;
01358 G4double ss;
01359
01360 while (itry1 < itry_max) {
01361 itry1++;
01362 ss = -std::log(inuclRndm());
01363 u = fmax * inuclRndm();
01364
01365 if (std::sqrt(ss) * std::exp(-ss) * (1.0 + ss/b) > u
01366 && ss < s4max) {
01367 ss = r0forAeq4 * std::sqrt(ss);
01368 coord1 = generateWithRandomAngles(ss).vect();
01369 coordinates.push_back(coord1);
01370
01371 if (verboseLevel > 2) {
01372 G4cout << " i " << i << " r " << coord1.mag() << G4endl;
01373 }
01374
01375 break;
01376 }
01377 }
01378
01379 if (itry1 == itry_max) {
01380 coord1.set(10000.,10000.,10000.);
01381 coordinates.push_back(coord1);
01382 break;
01383 }
01384 }
01385
01386 coord1 *= 0.0;
01387 for(G4int j = 0; j < ab -1; j++) coord1 -= coordinates[j];
01388
01389 coordinates.push_back(coord1);
01390
01391 if (verboseLevel > 2){
01392 G4cout << " last r " << coord1.mag() << G4endl;
01393 }
01394
01395 G4bool large_dist = false;
01396
01397 for (i = 0; i < ab-1; i++) {
01398 for (G4int j = i+1; j < ab; j++) {
01399
01400 G4double r2 = (coordinates[i]-coordinates[j]).mag2();
01401
01402 if (verboseLevel > 2){
01403 G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
01404 }
01405
01406 if (r2 > r_large2for4) {
01407 large_dist = true;
01408
01409 break;
01410 }
01411 }
01412
01413 if (large_dist) break;
01414 }
01415
01416 if (!large_dist) badco = false;
01417 }
01418 }
01419
01420 if(badco) {
01421 G4cout << " can not generate the nucleons coordinates for a "
01422 << ab << G4endl;
01423
01424 casparticles.clear();
01425 particles.clear();
01426 return;
01427
01428 } else {
01429 G4double p, u, x;
01430 G4LorentzVector mom;
01431
01432
01433 for (G4int i = 0; i < ab - 1; i++) {
01434 G4int itry2 = 0;
01435
01436 while(itry2 < itry_max) {
01437 itry2++;
01438 u = -std::log(0.879853 - 0.8798502 * inuclRndm());
01439 x = u * std::exp(-u);
01440
01441 if(x > inuclRndm()) {
01442 p = std::sqrt(0.01953 * u);
01443 mom = generateWithRandomAngles(p, massb);
01444 momentums.push_back(mom);
01445
01446 break;
01447 }
01448 }
01449
01450 if(itry2 == itry_max) {
01451 G4cout << " can not generate proper momentum for a "
01452 << ab << G4endl;
01453
01454 casparticles.clear();
01455 particles.clear();
01456 return;
01457 }
01458 }
01459
01460
01461 mom *= 0.;
01462 mom.setE(bullet->getEnergy()+target->getEnergy());
01463
01464 for(G4int j=0; j< ab-1; j++) mom -= momentums[j];
01465
01466 momentums.push_back(mom);
01467 }
01468 }
01469
01470
01471 G4double rb = 0.0;
01472 G4int i(0);
01473
01474 for(i = 0; i < G4int(coordinates.size()); i++) {
01475 G4double rp = coordinates[i].mag2();
01476
01477 if(rp > rb) rb = rp;
01478 }
01479
01480
01481 G4double s1 = std::sqrt(inuclRndm());
01482 G4double phi = randomPHI();
01483 G4double rz = (nuclei_radius + rb) * s1;
01484 G4ThreeVector global_pos(rz*std::cos(phi), rz*std::sin(phi),
01485 -(nuclei_radius+rb)*std::sqrt(1.0-s1*s1));
01486
01487 for (i = 0; i < G4int(coordinates.size()); i++) {
01488 coordinates[i] += global_pos;
01489 }
01490
01491
01492 raw_particles.clear();
01493
01494 for (G4int ipa = 0; ipa < ab; ipa++) {
01495 G4int knd = ipa < zb ? 1 : 2;
01496 raw_particles.push_back(G4InuclElementaryParticle(momentums[ipa], knd));
01497 }
01498
01499 G4InuclElementaryParticle dummy(small_ekin, 1);
01500 G4LorentzConvertor toTheBulletRestFrame(&dummy, bullet);
01501 toTheBulletRestFrame.toTheTargetRestFrame();
01502
01503 particleIterator ipart;
01504
01505 for (ipart = raw_particles.begin(); ipart != raw_particles.end(); ipart++) {
01506 ipart->setMomentum(toTheBulletRestFrame.backToTheLab(ipart->getMomentum()));
01507 }
01508
01509
01510
01511 for(G4int ip = 0; ip < G4int(raw_particles.size()); ip++) {
01512 G4LorentzVector mom = raw_particles[ip].getMomentum();
01513 G4double pmod = mom.rho();
01514 G4double t0 = -mom.vect().dot(coordinates[ip]) / pmod;
01515 G4double det = t0 * t0 + nuclei_radius * nuclei_radius
01516 - coordinates[ip].mag2();
01517 G4double tr = -1.0;
01518
01519 if(det > 0.0) {
01520 G4double t1 = t0 + std::sqrt(det);
01521 G4double t2 = t0 - std::sqrt(det);
01522
01523 if(std::fabs(t1) <= std::fabs(t2)) {
01524 if(t1 > 0.0) {
01525 if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
01526 }
01527
01528 if(tr < 0.0 && t2 > 0.0) {
01529
01530 if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
01531 }
01532
01533 } else {
01534 if(t2 > 0.0) {
01535
01536 if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
01537 }
01538
01539 if(tr < 0.0 && t1 > 0.0) {
01540 if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
01541 }
01542 }
01543
01544 }
01545
01546 if(tr >= 0.0) {
01547 coordinates[ip] += mom.vect()*tr / pmod;
01548 casparticles.push_back(G4CascadParticle(raw_particles[ip],
01549 coordinates[ip],
01550 number_of_zones, large, 0));
01551
01552 } else {
01553 particles.push_back(raw_particles[ip]);
01554 }
01555 }
01556 }
01557
01558 if(casparticles.size() == 0) {
01559 particles.clear();
01560
01561 G4cout << " can not generate proper distribution for " << itry_max
01562 << " steps " << G4endl;
01563 }
01564 }
01565 }
01566
01567 if(verboseLevel > 2){
01568 G4int ip(0);
01569
01570 G4cout << " cascad particles: " << casparticles.size() << G4endl;
01571 for(ip = 0; ip < G4int(casparticles.size()); ip++)
01572 G4cout << casparticles[ip] << G4endl;
01573
01574 G4cout << " outgoing particles: " << particles.size() << G4endl;
01575 for(ip = 0; ip < G4int(particles.size()); ip++)
01576 G4cout << particles[ip] << G4endl;
01577 }
01578
01579 return;
01580 }
01581
01582
01583
01584 G4double
01585 G4NucleiModel::inverseMeanFreePath(const G4CascadParticle& cparticle,
01586 const G4InuclElementaryParticle& target) {
01587 G4int ptype = cparticle.getParticle().type();
01588 G4int ip = target.type();
01589 G4int zone = cparticle.getCurrentZone();
01590
01591
01592 dummy_convertor.setBullet(cparticle.getParticle());
01593 dummy_convertor.setTarget(&target);
01594 G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
01595
01596
01597 G4double csec = (ip < 100) ? totalCrossSection(ekin, ptype*ip)
01598 : absorptionCrossSection(ekin, ptype);
01599
01600 if(verboseLevel > 2) {
01601 G4cout << " ip " << ip << " ekin " << ekin << " csec " << csec << G4endl;
01602 }
01603
01604 if (csec <= 0.) return 0.;
01605
01606
01607 return csec * getCurrentDensity(ip, zone);
01608 }
01609
01610
01611
01612 G4double
01613 G4NucleiModel::generateInteractionLength(const G4CascadParticle& cparticle,
01614 G4double path, G4double invmfp) const {
01615
01616 const G4double young_cut = std::sqrt(10.0) * 0.25;
01617 const G4double huge_num = 50.0;
01618 const G4double small = 1.0e-10;
01619
01620 if (invmfp < small) return path;
01621
01622 G4double pw = -path * invmfp;
01623 if (pw < -huge_num) pw = -huge_num;
01624 pw = 1.0 - std::exp(pw);
01625
01626 if (verboseLevel > 2)
01627 G4cout << " mfp " << 1./invmfp << " pw " << pw << G4endl;
01628
01629 G4double spath = path;
01630
01631
01632 if (isProjectile(cparticle) || (inuclRndm() < pw)) {
01633 spath = -std::log(1.0 - pw * inuclRndm()) / invmfp;
01634 if (cparticle.young(young_cut, spath)) spath = path;
01635
01636 if (verboseLevel > 2)
01637 G4cout << " spath " << spath << " path " << path << G4endl;
01638 }
01639
01640 return spath;
01641 }
01642
01643
01644
01645
01646 G4double G4NucleiModel::absorptionCrossSection(G4double ke, G4int type) const {
01647 if (type != pionPlus && type != pionMinus && type != pionZero &&
01648 type != photon) {
01649 G4cerr << " absorptionCrossSection only valid for incident pions" << G4endl;
01650 return 0.;
01651 }
01652
01653 G4double csec = 0.;
01654
01655
01656 if (type == pionPlus || type == pionMinus || type == pionZero) {
01657 if (ke < 0.3) csec = (0.1106 / std::sqrt(ke) - 0.8
01658 + 0.08 / ((ke-0.123)*(ke-0.123) + 0.0056) );
01659 else if (ke < 1.0) csec = 3.6735 * (1.0-ke)*(1.0-ke);
01660 }
01661
01662
01663 if (type == photon) {
01664 static const G4double gammaQDscale = G4CascadeParameters::gammaQDScale();
01665 static const G4double kebins[] =
01666 { 0.0, 0.01, 0.013, 0.018, 0.024, 0.032, 0.042, 0.056, 0.075, 0.1,
01667 0.13, 0.18, 0.24, 0.32, 0.42, 0.56, 0.75, 1.0, 1.3, 1.8,
01668 2.4, 3.2, 4.2, 5.6, 7.5, 10.0, 13.0, 18.0, 24.0, 32.0 };
01669 static const G4double gammaD[] =
01670 { 2.8, 1.3, 0.89, 0.56, 0.38, 0.27, 0.19, 0.14, 0.098,
01671 0.071, 0.054, 0.0003, 0.0007, 0.0027, 0.0014, 0.001, 0.0012, 0.0005,
01672 0.0003, 0.0002,0.0002, 0.0002, 0.0002, 0.0002, 0.0001, 0.0001, 0.0001,
01673 0.0001, 0.0001, 0.0001 };
01674 static G4CascadeInterpolator<30> interp(kebins);
01675 csec = interp.interpolate(ke, gammaD) * gammaQDscale;
01676 }
01677
01678 if (csec < 0.0) csec = 0.0;
01679
01680 if (verboseLevel > 2) {
01681 G4cout << " ekin " << ke << " abs. csec " << csec << " mb" << G4endl;
01682 }
01683
01684 return crossSectionUnits * csec;
01685 }
01686
01687
01688 G4double G4NucleiModel::totalCrossSection(G4double ke, G4int rtype) const
01689 {
01690
01691 const G4CascadeChannel* xsecTable = G4CascadeChannelTables::GetTable(rtype);
01692 if (!xsecTable) {
01693 G4cerr << " unknown collison type = " << rtype << G4endl;
01694 return 0.;
01695 }
01696
01697 return (crossSectionUnits * xsecTable->getCrossSection(ke));
01698 }