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 #include <algorithm>
00115
00116 #include "G4IntraNucleiCascader.hh"
00117 #include "G4SystemOfUnits.hh"
00118 #include "G4CascadeChannelTables.hh"
00119 #include "G4CascadeCoalescence.hh"
00120 #include "G4CascadeParameters.hh"
00121 #include "G4CascadeRecoilMaker.hh"
00122 #include "G4CascadParticle.hh"
00123 #include "G4CollisionOutput.hh"
00124 #include "G4DecayProducts.hh"
00125 #include "G4DecayTable.hh"
00126 #include "G4ElementaryParticleCollider.hh"
00127 #include "G4ExitonConfiguration.hh"
00128 #include "G4HadTmpUtil.hh"
00129 #include "G4InuclElementaryParticle.hh"
00130 #include "G4InuclNuclei.hh"
00131 #include "G4InuclParticleNames.hh"
00132 #include "G4InuclSpecialFunctions.hh"
00133 #include "G4KineticTrack.hh"
00134 #include "G4KineticTrackVector.hh"
00135 #include "G4LorentzConvertor.hh"
00136 #include "G4Neutron.hh"
00137 #include "G4NucleiModel.hh"
00138 #include "G4ParticleLargerEkin.hh"
00139 #include "G4Proton.hh"
00140 #include "G4V3DNucleus.hh"
00141 #include "Randomize.hh"
00142
00143 using namespace G4InuclParticleNames;
00144 using namespace G4InuclSpecialFunctions;
00145
00146
00147
00148 const G4int G4IntraNucleiCascader::itry_max = 100;
00149 const G4int G4IntraNucleiCascader::reflection_cut = 50;
00150 const G4double G4IntraNucleiCascader::small_ekin = 0.001*MeV;
00151 const G4double G4IntraNucleiCascader::quasielast_cut = 1*MeV;
00152
00153
00154 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
00155
00156 G4IntraNucleiCascader::G4IntraNucleiCascader()
00157 : G4CascadeColliderBase("G4IntraNucleiCascader"), model(new G4NucleiModel),
00158 theElementaryParticleCollider(new G4ElementaryParticleCollider),
00159 theRecoilMaker(new G4CascadeRecoilMaker), theClusterMaker(0),
00160 tnuclei(0), bnuclei(0), bparticle(0),
00161 minimum_recoil_A(0.), coulombBarrier(0.),
00162 nucleusTarget(new G4InuclNuclei),
00163 protonTarget(new G4InuclElementaryParticle) {
00164 if (G4CascadeParameters::doCoalescence())
00165 theClusterMaker = new G4CascadeCoalescence;
00166 }
00167
00168 G4IntraNucleiCascader::~G4IntraNucleiCascader() {
00169 delete model;
00170 delete theElementaryParticleCollider;
00171 delete theRecoilMaker;
00172 delete theClusterMaker;
00173 delete nucleusTarget;
00174 delete protonTarget;
00175 }
00176
00177 void G4IntraNucleiCascader::setVerboseLevel(G4int verbose) {
00178 G4CascadeColliderBase::setVerboseLevel(verbose);
00179 model->setVerboseLevel(verbose);
00180 theElementaryParticleCollider->setVerboseLevel(verbose);
00181 theRecoilMaker->setVerboseLevel(verbose);
00182 }
00183
00184
00185
00186 void G4IntraNucleiCascader::collide(G4InuclParticle* bullet,
00187 G4InuclParticle* target,
00188 G4CollisionOutput& globalOutput) {
00189 if (verboseLevel) G4cout << " >>> G4IntraNucleiCascader::collide " << G4endl;
00190
00191 if (!initialize(bullet, target)) return;
00192
00193 G4int itry = 0;
00194 do {
00195 newCascade(++itry);
00196 setupCascade();
00197 generateCascade();
00198 } while (!finishCascade() && itry<itry_max);
00199
00200 finalize(itry, bullet, target, globalOutput);
00201 }
00202
00203
00204
00205
00206 void G4IntraNucleiCascader::rescatter(G4InuclParticle* bullet,
00207 G4KineticTrackVector* theSecondaries,
00208 G4V3DNucleus* theNucleus,
00209 G4CollisionOutput& globalOutput) {
00210 if (verboseLevel)
00211 G4cout << " >>> G4IntraNucleiCascader::rescatter " << G4endl;
00212
00213 G4InuclParticle* target = createTarget(theNucleus);
00214 if (!initialize(bullet, target)) return;
00215
00216 G4int itry = 0;
00217 do {
00218 newCascade(++itry);
00219 preloadCascade(theNucleus, theSecondaries);
00220 generateCascade();
00221 } while (!finishCascade() && itry<itry_max);
00222
00223 finalize(itry, bullet, target, globalOutput);
00224 }
00225
00226
00227 G4bool G4IntraNucleiCascader::initialize(G4InuclParticle* bullet,
00228 G4InuclParticle* target) {
00229 if (verboseLevel>1)
00230 G4cout << " >>> G4IntraNucleiCascader::initialize " << G4endl;
00231
00232
00233 theRecoilMaker->setTolerance(small_ekin);
00234
00235 interCase.set(bullet,target);
00236
00237 if (verboseLevel > 3) {
00238 G4cout << *interCase.getBullet() << G4endl
00239 << *interCase.getTarget() << G4endl;
00240 }
00241
00242
00243 bnuclei = dynamic_cast<G4InuclNuclei*>(interCase.getBullet());
00244 bparticle = dynamic_cast<G4InuclElementaryParticle*>(interCase.getBullet());
00245
00246 if (!bnuclei && !bparticle) {
00247 G4cerr << " G4IntraNucleiCascader: projectile is not a valid particle."
00248 << G4endl;
00249 return false;
00250 }
00251
00252
00253 tnuclei = dynamic_cast<G4InuclNuclei*>(interCase.getTarget());
00254 if (!tnuclei) {
00255 if (verboseLevel)
00256 G4cerr << " Target is not a nucleus. Abandoning." << G4endl;
00257 return false;
00258 }
00259
00260 model->generateModel(tnuclei);
00261 coulombBarrier = 0.00126*tnuclei->getZ() / (1.+G4cbrt(tnuclei->getA()));
00262
00263
00264
00265 minimum_recoil_A = 0.;
00266
00267 if (verboseLevel > 3) {
00268 G4LorentzVector momentum_in = bullet->getMomentum() + target->getMomentum();
00269 G4cout << " intitial momentum E " << momentum_in.e() << " Px "
00270 << momentum_in.x() << " Py " << momentum_in.y() << " Pz "
00271 << momentum_in.z() << G4endl;
00272 }
00273
00274 return true;
00275 }
00276
00277
00278
00279 void G4IntraNucleiCascader::newCascade(G4int itry) {
00280 if (verboseLevel > 1) {
00281 G4cout << " IntraNucleiCascader itry " << itry << " inter_case "
00282 << interCase.code() << G4endl;
00283 }
00284
00285 model->reset();
00286 output.reset();
00287 new_cascad_particles.clear();
00288 theExitonConfiguration.clear();
00289
00290 cascad_particles.clear();
00291 }
00292
00293
00294
00295
00296 void G4IntraNucleiCascader::setupCascade() {
00297 if (verboseLevel > 1)
00298 G4cout << " >>> G4IntraNucleiCascader::setupCascade" << G4endl;
00299
00300 if (interCase.hadNucleus()) {
00301 if (verboseLevel > 3)
00302 G4cout << " bparticle charge " << bparticle->getCharge()
00303 << " baryon number " << bparticle->baryon() << G4endl;
00304
00305 cascad_particles.push_back(model->initializeCascad(bparticle));
00306 } else {
00307 G4int ab = bnuclei->getA();
00308 G4int zb = bnuclei->getZ();
00309
00310 G4NucleiModel::modelLists all_particles;
00311 model->initializeCascad(bnuclei, tnuclei, all_particles);
00312
00313 cascad_particles = all_particles.first;
00314 output.addOutgoingParticles(all_particles.second);
00315
00316 if (cascad_particles.size() == 0) {
00317 G4int i;
00318
00319 for (i = 0; i < ab; i++) {
00320 G4int knd = i < zb ? 1 : 2;
00321 theExitonConfiguration.incrementQP(knd);
00322 };
00323
00324 G4int ihn = G4int(2 * (ab-zb) * inuclRndm() + 0.5);
00325 G4int ihz = G4int(2 * zb * inuclRndm() + 0.5);
00326
00327 for (i = 0; i < ihn; i++) theExitonConfiguration.incrementHoles(2);
00328 for (i = 0; i < ihz; i++) theExitonConfiguration.incrementHoles(1);
00329 }
00330 }
00331 }
00332
00333
00334
00335
00336 void G4IntraNucleiCascader::generateCascade() {
00337 if (verboseLevel>1) G4cout << " generateCascade " << G4endl;
00338
00339 G4int iloop = 0;
00340 while (!cascad_particles.empty() && !model->empty()) {
00341 iloop++;
00342
00343 if (verboseLevel > 2) {
00344 G4cout << " Iteration " << iloop << ": Number of cparticles "
00345 << cascad_particles.size() << " last one: \n"
00346 << cascad_particles.back() << G4endl;
00347 }
00348
00349 model->generateParticleFate(cascad_particles.back(),
00350 theElementaryParticleCollider,
00351 new_cascad_particles);
00352
00353 if (verboseLevel > 2) {
00354 G4cout << " After generate fate: New particles "
00355 << new_cascad_particles.size() << G4endl
00356 << " Discarding last cparticle from list " << G4endl;
00357 }
00358
00359 cascad_particles.pop_back();
00360
00361
00362
00363 if (new_cascad_particles.size() == 1) {
00364 const G4CascadParticle& currentCParticle = new_cascad_particles[0];
00365
00366 if (model->stillInside(currentCParticle)) {
00367 if (verboseLevel > 3)
00368 G4cout << " particle still inside nucleus " << G4endl;
00369
00370 if (currentCParticle.getNumberOfReflections() < reflection_cut &&
00371 model->worthToPropagate(currentCParticle)) {
00372 if (verboseLevel > 3) G4cout << " continue reflections " << G4endl;
00373 cascad_particles.push_back(currentCParticle);
00374 } else {
00375 processTrappedParticle(currentCParticle);
00376 }
00377
00378 } else {
00379 if (verboseLevel > 3) G4cout << " possible escape " << G4endl;
00380
00381 const G4InuclElementaryParticle& currentParticle =
00382 currentCParticle.getParticle();
00383
00384 G4double KE = currentParticle.getKineticEnergy();
00385 G4double mass = currentParticle.getMass();
00386 G4double Q = currentParticle.getCharge();
00387
00388 if (verboseLevel > 3)
00389 G4cout << " KE " << KE << " barrier " << Q*coulombBarrier << G4endl;
00390
00391 if (KE < Q*coulombBarrier) {
00392
00393 G4double CBP = 0.0;
00394
00395
00396
00397 if (KE > 0.0001) CBP = std::exp(-0.0181*0.5*tnuclei->getZ()*
00398 (1./KE - 1./coulombBarrier)*
00399 std::sqrt(mass*(coulombBarrier-KE)) );
00400
00401 if (G4UniformRand() < CBP) {
00402 if (verboseLevel > 3)
00403 G4cout << " tunneled\n" << currentParticle << G4endl;
00404
00405
00406 output.addOutgoingParticle(currentParticle);
00407 } else {
00408 processTrappedParticle(currentCParticle);
00409 }
00410 } else {
00411 output.addOutgoingParticle(currentParticle);
00412
00413 if (verboseLevel > 3)
00414 G4cout << " Goes out\n" << output.getOutgoingParticles().back()
00415 << G4endl;
00416 }
00417 }
00418 } else {
00419 if (verboseLevel > 3)
00420 G4cout << " interacted, adding new to list " << G4endl;
00421
00422 cascad_particles.insert(cascad_particles.end(),
00423 new_cascad_particles.begin(),
00424 new_cascad_particles.end());
00425
00426 std::pair<G4int, G4int> holes = model->getTypesOfNucleonsInvolved();
00427 if (verboseLevel > 3)
00428 G4cout << " adding new exciton holes " << holes.first << ","
00429 << holes.second << G4endl;
00430
00431 theExitonConfiguration.incrementHoles(holes.first);
00432
00433 if (holes.second > 0)
00434 theExitonConfiguration.incrementHoles(holes.second);
00435 }
00436
00437
00438 theRecoilMaker->collide(interCase.getBullet(), interCase.getTarget(),
00439 output, cascad_particles);
00440
00441 G4double aresid = theRecoilMaker->getRecoilA();
00442 if (verboseLevel > 2) {
00443 G4cout << " cparticles remaining " << cascad_particles.size()
00444 << " nucleus (model) has "
00445 << model->getNumberOfNeutrons() << " n, "
00446 << model->getNumberOfProtons() << " p "
00447 << " residual fragment A " << aresid << G4endl;
00448 }
00449
00450 if (aresid <= minimum_recoil_A) return;
00451 }
00452 }
00453
00454
00455
00456
00457 G4bool G4IntraNucleiCascader::finishCascade() {
00458 if (verboseLevel > 1)
00459 G4cout << " >>> G4IntraNucleiCascader::finishCascade ?" << G4endl;
00460
00461
00462 output.addOutgoingParticles(cascad_particles);
00463 cascad_particles.clear();
00464
00465
00466 if (verboseLevel>3) {
00467 G4cout << " G4IntraNucleiCascader finished" << G4endl;
00468 output.printCollisionOutput();
00469 }
00470
00471
00472 if (theClusterMaker) {
00473 theClusterMaker->setVerboseLevel(verboseLevel);
00474 theClusterMaker->FindClusters(output);
00475
00476
00477 if (verboseLevel>3) G4cout << " Recomputing recoil fragment" << G4endl;
00478 theRecoilMaker->collide(interCase.getBullet(), interCase.getTarget(),
00479 output);
00480 if (verboseLevel>3) {
00481 G4cout << " After cluster coalescence" << G4endl;
00482 output.printCollisionOutput();
00483 }
00484 }
00485
00486
00487 G4int afin = theRecoilMaker->getRecoilA();
00488 G4int zfin = theRecoilMaker->getRecoilZ();
00489
00490
00491
00492
00493 if (!theRecoilMaker->goodFragment() && !theRecoilMaker->wholeEvent()) {
00494 if (verboseLevel > 1)
00495 G4cerr << " Recoil nucleus is not physical: A=" << afin << " Z="
00496 << zfin << G4endl;
00497 return false;
00498 }
00499
00500 const G4LorentzVector& presid = theRecoilMaker->getRecoilMomentum();
00501
00502 if (verboseLevel > 1) {
00503 G4cout << " afin " << afin << " zfin " << zfin << G4endl;
00504 }
00505
00506 if (afin == 0) return true;
00507
00508 if (afin == 1) {
00509 G4int last_type = (zfin==1) ? 1 : 2;
00510
00511 G4double mass = G4InuclElementaryParticle::getParticleMass(last_type);
00512 G4double mres = presid.m();
00513
00514
00515 if (mres-mass < -small_ekin) {
00516 if (verboseLevel > 2) G4cerr << " unphysical recoil nucleon" << G4endl;
00517 return false;
00518 }
00519
00520 if (mres-mass > small_ekin) {
00521 if (verboseLevel > 2)
00522 G4cerr << " extra energy with recoil nucleon" << G4endl;
00523
00524
00525
00526 }
00527
00528 G4InuclElementaryParticle last_particle(presid, last_type,
00529 G4InuclParticle::INCascader);
00530
00531 if (verboseLevel > 3) {
00532 G4cout << " adding recoiling nucleon to output list\n"
00533 << last_particle << G4endl;
00534 }
00535
00536 output.addOutgoingParticle(last_particle);
00537
00538
00539 theRecoilMaker->collide(interCase.getBullet(), interCase.getTarget(),
00540 output);
00541 }
00542
00543
00544 if (output.numberOfOutgoingParticles() == 1) {
00545 G4double Eex = theRecoilMaker->getRecoilExcitation();
00546 if (std::abs(Eex) < quasielast_cut) {
00547 if (verboseLevel > 3) {
00548 G4cout << " quasi-elastic scatter with " << Eex << " MeV recoil"
00549 << G4endl;
00550 }
00551
00552 theRecoilMaker->setRecoilExcitation(Eex=0.);
00553 if (verboseLevel > 3) {
00554 G4cout << " Eex reset to " << theRecoilMaker->getRecoilExcitation()
00555 << G4endl;
00556 }
00557 }
00558 }
00559
00560 if (theRecoilMaker->goodNucleus()) {
00561 theRecoilMaker->addExcitonConfiguration(theExitonConfiguration);
00562
00563 G4Fragment* recoilFrag = theRecoilMaker->makeRecoilFragment();
00564 if (!recoilFrag) {
00565 G4cerr << "Got null pointer for recoil fragment!" << G4endl;
00566 return false;
00567 }
00568
00569 if (verboseLevel > 2)
00570 G4cout << " adding recoil fragment to output list" << G4endl;
00571
00572 output.addRecoilFragment(*recoilFrag);
00573 }
00574
00575
00576 std::vector<G4InuclElementaryParticle>& opart = output.getOutgoingParticles();
00577 std::sort(opart.begin(), opart.end(), G4ParticleLargerEkin());
00578
00579
00580 if (theRecoilMaker->wholeEvent() || theRecoilMaker->goodNucleus()) {
00581 output.setVerboseLevel(verboseLevel);
00582 output.setOnShell(interCase.getBullet(), interCase.getTarget());
00583 output.setVerboseLevel(0);
00584
00585 if (output.acceptable()) return true;
00586 else if (verboseLevel>2) G4cerr << " Cascade setOnShell failed." << G4endl;
00587 }
00588
00589
00590 if (afin <= minimum_recoil_A && minimum_recoil_A < tnuclei->getA()) {
00591 ++minimum_recoil_A;
00592 if (verboseLevel > 3) {
00593 G4cout << " minimum recoil fragment increased to A " << minimum_recoil_A
00594 << G4endl;
00595 }
00596 }
00597
00598 if (verboseLevel>2) G4cerr << " Cascade failed. Retrying..." << G4endl;
00599 return false;
00600 }
00601
00602
00603
00604
00605 void
00606 G4IntraNucleiCascader::finalize(G4int itry, G4InuclParticle* bullet,
00607 G4InuclParticle* target,
00608 G4CollisionOutput& globalOutput) {
00609 if (itry >= itry_max) {
00610 if (verboseLevel) {
00611 G4cout << " IntraNucleiCascader-> no inelastic interaction after "
00612 << itry << " attempts " << G4endl;
00613 }
00614
00615 output.trivialise(bullet, target);
00616 } else if (verboseLevel) {
00617 G4cout << " IntraNucleiCascader output after trials " << itry << G4endl;
00618 }
00619
00620
00621 globalOutput.add(output);
00622 }
00623
00624
00625
00626
00627 G4InuclParticle*
00628 G4IntraNucleiCascader::createTarget(G4V3DNucleus* theNucleus) {
00629 G4int theNucleusA = theNucleus->GetMassNumber();
00630 G4int theNucleusZ = theNucleus->GetCharge();
00631
00632 if (theNucleusA > 1) {
00633 if (!nucleusTarget) nucleusTarget = new G4InuclNuclei;
00634 nucleusTarget->fill(0., theNucleusA, theNucleusZ, 0.);
00635 return nucleusTarget;
00636 } else {
00637 if (!protonTarget) protonTarget = new G4InuclElementaryParticle;
00638 protonTarget->fill(0., (theNucleusZ==1)?proton:neutron);
00639 return protonTarget;
00640 }
00641
00642 return 0;
00643 }
00644
00645
00646
00647 void
00648 G4IntraNucleiCascader::preloadCascade(G4V3DNucleus* theNucleus,
00649 G4KineticTrackVector* theSecondaries) {
00650 if (verboseLevel > 1)
00651 G4cout << " >>> G4IntraNucleiCascader::preloadCascade" << G4endl;
00652
00653 copyWoundedNucleus(theNucleus);
00654 copySecondaries(theSecondaries);
00655 }
00656
00657 void G4IntraNucleiCascader::copyWoundedNucleus(G4V3DNucleus* theNucleus) {
00658 if (verboseLevel > 1)
00659 G4cout << " >>> G4IntraNucleiCascader::copyWoundedNucleus" << G4endl;
00660
00661
00662 theExitonConfiguration.clear();
00663 hitNucleons.clear();
00664 if (theNucleus->StartLoop()) {
00665 G4Nucleon* nucl = 0;
00666 G4int nuclType = 0;
00667 while ((nucl = theNucleus->GetNextNucleon())) {
00668 if (nucl->AreYouHit()) {
00669 nuclType = G4InuclElementaryParticle::type(nucl->GetParticleType());
00670 theExitonConfiguration.incrementHoles(nuclType);
00671 hitNucleons.push_back(nucl->GetPosition());
00672 }
00673 }
00674 }
00675
00676 if (verboseLevel > 3)
00677 G4cout << " nucleus has " << theExitonConfiguration.neutronHoles
00678 << " neutrons hit, " << theExitonConfiguration.protonHoles
00679 << " protons hit" << G4endl;
00680
00681
00682 model->reset(theExitonConfiguration.neutronHoles,
00683 theExitonConfiguration.protonHoles, &hitNucleons);
00684 }
00685
00686 void
00687 G4IntraNucleiCascader::copySecondaries(G4KineticTrackVector* secondaries) {
00688 if (verboseLevel > 1)
00689 G4cout << " >>> G4IntraNucleiCascader::copySecondaries" << G4endl;
00690
00691 for (size_t i=0; i<secondaries->size(); i++) {
00692 if (verboseLevel > 3) G4cout << " processing secondary " << i << G4endl;
00693
00694 processSecondary((*secondaries)[i]);
00695 }
00696
00697
00698 std::sort(cascad_particles.begin(), cascad_particles.end(),
00699 G4ParticleLargerEkin());
00700
00701 if (verboseLevel > 2) {
00702 G4cout << " Original list of " << secondaries->size() << " secondaries"
00703 << " produced " << cascad_particles.size() << " cascade, "
00704 << output.numberOfOutgoingParticles() << " released particles, "
00705 << output.numberOfOutgoingNuclei() << " fragments" << G4endl;
00706 }
00707 }
00708
00709
00710
00711
00712 void G4IntraNucleiCascader::processSecondary(const G4KineticTrack* ktrack) {
00713 if (!ktrack) return;
00714
00715
00716 G4ParticleDefinition* kpd = ktrack->GetDefinition();
00717 if (!kpd) return;
00718
00719 G4int ktype = G4InuclElementaryParticle::type(kpd);
00720 if (!ktype) {
00721 releaseSecondary(ktrack);
00722 return;
00723 }
00724
00725 if (verboseLevel > 1) {
00726 G4cout << " >>> G4IntraNucleiCascader::processSecondary "
00727 << kpd->GetParticleName() << G4endl;
00728 }
00729
00730
00731 cascad_particles.resize(cascad_particles.size()+1);
00732 G4CascadParticle& cpart = cascad_particles.back();
00733
00734
00735 cpart.getParticle().fill(ktrack->Get4Momentum()/GeV, ktype);
00736 cpart.setGeneration(0);
00737 cpart.setMovingInsideNuclei();
00738 cpart.initializePath(0);
00739
00740
00741 G4ThreeVector cpos = ktrack->GetPosition()/model->getRadiusUnits();
00742
00743 cpart.updatePosition(cpos);
00744 cpart.updateZone(model->getZone(cpos.mag()));
00745
00746 if (verboseLevel > 2)
00747 G4cout << " Created cascade particle \n" << cpart << G4endl;
00748 }
00749
00750
00751
00752
00753 void G4IntraNucleiCascader::releaseSecondary(const G4KineticTrack* ktrack) {
00754 G4ParticleDefinition* kpd = ktrack->GetDefinition();
00755
00756 if (verboseLevel > 1) {
00757 G4cout << " >>> G4IntraNucleiCascader::releaseSecondary "
00758 << kpd->GetParticleName() << G4endl;
00759 }
00760
00761
00762 if (dynamic_cast<G4Ions*>(kpd)) {
00763
00764 output.getOutgoingNuclei().resize(output.numberOfOutgoingNuclei()+1);
00765 G4InuclNuclei& inucl = output.getOutgoingNuclei().back();
00766
00767 inucl.fill(ktrack->Get4Momentum()/GeV,
00768 kpd->GetAtomicMass(), kpd->GetAtomicNumber());
00769 if (verboseLevel > 2)
00770 G4cout << " Created pre-cascade fragment\n" << inucl << G4endl;
00771 } else {
00772
00773 output.getOutgoingParticles().resize(output.numberOfOutgoingParticles()+1);
00774 G4InuclElementaryParticle& ipart = output.getOutgoingParticles().back();
00775
00776
00777 ipart.fill(ktrack->Get4Momentum()/GeV, ktrack->GetDefinition());
00778 if (verboseLevel > 2)
00779 G4cout << " Created invalid pre-cascade particle\n" << ipart << G4endl;
00780 }
00781 }
00782
00783
00784
00785
00786 void G4IntraNucleiCascader::
00787 processTrappedParticle(const G4CascadParticle& trapped) {
00788 const G4InuclElementaryParticle& trappedP = trapped.getParticle();
00789
00790 G4int xtype = trappedP.type();
00791 if (verboseLevel > 3) G4cout << " exciton of type " << xtype << G4endl;
00792
00793 if (trappedP.nucleon()) {
00794 theExitonConfiguration.incrementQP(xtype);
00795 return;
00796 }
00797
00798 if (trappedP.hyperon()) {
00799 decayTrappedParticle(trapped);
00800 return;
00801 }
00802
00803
00804
00805 if (verboseLevel > 3) {
00806 G4cout << " non-standard should be absorbed, now released\n"
00807 << trapped << G4endl;
00808 }
00809
00810 output.addOutgoingParticle(trappedP);
00811 }
00812
00813
00814
00815
00816 void G4IntraNucleiCascader::
00817 decayTrappedParticle(const G4CascadParticle& trapped) {
00818 if (verboseLevel > 3)
00819 G4cout << " unstable must be decayed in flight" << G4endl;
00820
00821 const G4InuclElementaryParticle& trappedP = trapped.getParticle();
00822
00823 G4DecayTable* unstable = trappedP.getDefinition()->GetDecayTable();
00824 if (!unstable) {
00825 if (verboseLevel > 3)
00826 G4cerr << " no decay table! Releasing trapped particle" << G4endl;
00827
00828 output.addOutgoingParticle(trappedP);
00829 return;
00830 }
00831
00832
00833 G4DecayProducts* daughters = unstable->SelectADecayChannel()->DecayIt( trappedP.getDefinition()->GetPDGMass() );
00834 if (!daughters) {
00835 if (verboseLevel > 3)
00836 G4cerr << " no daughters! Releasing trapped particle" << G4endl;
00837
00838 output.addOutgoingParticle(trappedP);
00839 return;
00840 }
00841
00842 if (verboseLevel > 3)
00843 G4cout << " " << daughters->entries() << " decay daughters" << G4endl;
00844
00845
00846 G4double decayEnergy = trappedP.getEnergy();
00847 G4ThreeVector decayDir = trappedP.getMomentum().vect().unit();
00848 daughters->Boost(decayEnergy, decayDir);
00849
00850
00851 const G4ThreeVector& decayPos = trapped.getPosition();
00852 G4int zone = trapped.getCurrentZone();
00853 G4int gen = trapped.getGeneration()+1;
00854
00855 for (G4int i=0; i<daughters->entries(); i++) {
00856 G4DynamicParticle* idaug = (*daughters)[i];
00857
00858 G4InuclElementaryParticle idaugEP(*idaug, G4InuclParticle::INCascader);
00859
00860
00861 if (G4CascadeChannelTables::GetTable(idaugEP.type())) {
00862 if (verboseLevel > 3) G4cout << " propagating " << idaugEP << G4endl;
00863 cascad_particles.push_back(G4CascadParticle(idaugEP,decayPos,zone,0.,gen));
00864 } else {
00865 if (verboseLevel > 3) G4cout << " releasing " << idaugEP << G4endl;
00866 output.addOutgoingParticle(idaugEP);
00867 }
00868 }
00869 }