124 using namespace G4InuclParticleNames;
125 using namespace G4InuclSpecialFunctions;
140 G4cout <<
" >>> G4ElementaryParticleCollider::collide" <<
G4endl;
143 G4cerr <<
" ElementaryParticleCollider -> can collide only particle with particle "
148 #ifdef G4CASCADE_DEBUG_SAMPLER
149 static G4bool doPrintTables =
true;
152 doPrintTables =
false;
165 if (!particle1 || !particle2) {
166 G4cerr <<
" ElementaryParticleCollider -> can only collide hadrons"
176 G4cerr <<
" ElementaryParticleCollider -> cannot collide "
199 generateSCMfinalState(ekin, etot_scm, particle1, particle2);
206 G4cerr <<
" ElementaryParticleCollider -> can only collide pi,mu,gamma with"
207 <<
" dibaryons " <<
G4endl;
214 generateSCMmuonAbsorption(etot_scm, particle1, particle2);
216 generateSCMpionAbsorption(etot_scm, particle1, particle2);
220 if (particles.empty()) {
222 G4cerr <<
" ElementaryParticleCollider -> failed to collide "
233 for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
235 ipart->setMomentum(mom);
240 G4cout <<
" incoming particles: \n" << *particle1 << G4endl
241 << *particle2 << G4endl
242 <<
" outgoing particles: " <<
G4endl;
243 for(ipart = particles.begin(); ipart != particles.end(); ipart++)
244 G4cout << *ipart << G4endl;
246 G4cout <<
" <<< Non-conservation in G4ElementaryParticleCollider"
256 G4ElementaryParticleCollider::generateMultiplicity(
G4int is,
265 G4cerr <<
" G4ElementaryParticleCollider: Unknown interaction channel "
266 << is <<
" - multiplicity not generated " <<
G4endl;
270 G4cout <<
" G4ElementaryParticleCollider::generateMultiplicity: "
271 <<
" multiplicity = " << mul <<
G4endl;
279 G4ElementaryParticleCollider::generateOutgoingPartTypes(
G4int is,
G4int mult,
282 particle_kinds.clear();
289 G4cerr <<
" G4ElementaryParticleCollider: Unknown interaction channel "
290 << is <<
" - outgoing kinds not generated " <<
G4endl;
298 G4ElementaryParticleCollider::generateSCMfinalState(
G4double ekin,
303 G4cout <<
" >>> G4ElementaryParticleCollider::generateSCMfinalState"
309 const G4int itry_max = 10;
314 G4int is = type1 * type2;
318 G4int multiplicity = 0;
322 while (generate && itry++ < itry_max) {
324 particle_kinds.clear();
327 multiplicity = generateMultiplicity(is, ekin);
329 generateOutgoingPartTypes(is, multiplicity, ekin);
330 if (particle_kinds.empty()) {
332 G4cout <<
" generateOutgoingPartTypes failed mult " << multiplicity
338 fillOutgoingMasses();
341 fsGenerator.
Configure(particle1, particle2, particle_kinds);
342 generate = !fsGenerator.
Generate(etot_scm, masses, scm_momentums);
345 if (itry >= itry_max) {
347 G4cout <<
" generateSCMfinalState failed " << itry <<
" attempts"
354 particles.resize(multiplicity);
355 for (
G4int i=0; i<multiplicity; i++) {
356 particles[i].fill(scm_momentums[i], particle_kinds[i],
361 G4cout <<
" <<< G4ElementaryParticleCollider::generateSCMfinalState"
371 void G4ElementaryParticleCollider::fillOutgoingMasses() {
372 G4int mult = particle_kinds.size();
374 masses.resize(mult,0.);
375 masses2.resize(mult,0.);
377 for (
G4int i = 0; i < mult; i++) {
379 masses2[i] = masses[i] * masses[i];
388 G4ElementaryParticleCollider::generateSCMpionAbsorption(
G4double etot_scm,
392 G4cout <<
" >>> G4ElementaryParticleCollider::generateSCMpionAbsorption"
398 particle_kinds.clear();
405 G4cerr <<
" pion absorption: "
412 if (!splitQuasiDeuteron(type2))
return;
414 fillOutgoingMasses();
416 G4double a = 0.5 * (etot_scm * etot_scm - masses2[0] - masses2[1]);
418 G4double pmod = std::sqrt((a*a - masses2[0]*masses2[1])
419 / (masses2[0] + masses2[1] + 2.0*a));
433 G4ElementaryParticleCollider::generateSCMmuonAbsorption(
G4double etot_scm,
438 G4cout <<
" >>> G4ElementaryParticleCollider::generateSCMmuonAbsorption"
446 scm_momentums.clear();
447 scm_momentums.resize(3);
449 particle_kinds.clear();
458 G4cerr <<
" mu- absorption: "
465 if (!splitQuasiDeuteron(type2))
return;
466 particle_kinds.push_back(
mnu);
468 fillOutgoingMasses();
470 G4GDecay3 breakup(etot_scm, masses[0], masses[1], masses[2]);
471 std::vector<G4ThreeVector> theMomenta = breakup.GetThreeBodyMomenta();
473 for (
size_t i=0; i<3; i++) {
474 scm_momentums[i].setVectM(theMomenta[i], masses[i]);
482 G4bool G4ElementaryParticleCollider::splitQuasiDeuteron(
G4int qdtype) {
484 G4cerr <<
" type " << qdtype <<
" not dibaryon!" <<
G4endl;
488 G4int b2 = qdtype % 10;
489 G4int b1 = (qdtype/10) % 10;
491 particle_kinds.push_back(b1);
492 particle_kinds.push_back(b2);
static const G4CascadeChannel * GetTable(G4int initialState)
G4bool Generate(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
virtual G4bool useEPCollider(G4InuclParticle *bullet, G4InuclParticle *target) const
virtual void getOutgoingParticleTypes(std::vector< G4int > &kinds, G4int mult, G4double ke) const =0
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
static G4bool useQuasiDeuteron(G4int ptype, G4int qdtype=0)
static G4double getParticleMass(G4int type)
void setBullet(const G4InuclParticle *bullet)
G4double getTotalSCMEnergy() const
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
const G4String & GetParticleName() const
void setVectM(const Hep3Vector &spatial, double mass)
G4ElementaryParticleCollider()
G4GLOB_DLL std::ostream G4cout
void setVerbose(G4int vb=0)
void Configure(G4InuclElementaryParticle *bullet, G4InuclElementaryParticle *target, const std::vector< G4int > &particle_kinds)
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
G4double getKinEnergyInTheTRS() const
G4InteractionCase interCase
static void Print(std::ostream &os=G4cout)
void addOutgoingParticles(const std::vector< G4InuclElementaryParticle > &particles)
void set(G4InuclParticle *part1, G4InuclParticle *part2)
void SetVerboseLevel(G4int verbose)
virtual G4bool validateOutput(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
std::vector< G4InuclElementaryParticle >::iterator particleIterator
G4bool quasi_deutron() const
virtual G4int getMultiplicity(G4double ke) const =0
G4ParticleDefinition * getDefinition() const
G4double getMomModule() const
G4GLOB_DLL std::ostream G4cerr
G4bool isNeutrino() const
void setTarget(const G4InuclParticle *target)