Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ElementaryParticleCollider.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4ElementaryParticleCollider.cc 71940 2013-06-28 19:04:44Z mkelsey $
27 //
28 // 20100114 M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
29 // 20100316 D. Wright (restored by M. Kelsey) -- Replace original (incorrect)
30 // pp, pn, nn 2-body to 2-body scattering angular distributions
31 // with a new parametrization of elastic scattering data using
32 // the sum of two exponentials.
33 // 20100319 M. Kelsey -- Use new generateWithRandomAngles for theta,phi stuff;
34 // eliminate some unnecessary std::pow()
35 // 20100407 M. Kelsey -- Replace std::vector<>::resize(0) with ::clear()
36 // Eliminate return-by-value std::vector<> by creating buffers
37 // buffers for particles, momenta, and particle types.
38 // The following functions now return void and are non-const:
39 // ::generateSCMfinalState()
40 // ::generateMomModules()
41 // ::generateStrangeChannelPartTypes()
42 // ::generateSCMpionAbsorption()
43 // 20100408 M. Kelsey -- Follow changes to G4*Sampler to pass particle_kinds
44 // as input buffer.
45 // 20100413 M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
46 // 20100428 M. Kelsey -- Use G4InuclParticleNames enum
47 // 20100429 M. Kelsey -- Change "photon()" to "isPhoton()"
48 // 20100507 M. Kelsey -- Rationalize multiplicity returns to be actual value
49 // 20100511 M. Kelsey -- Replace G4PionSampler and G4NucleonSampler with new
50 // pi-N and N-N classes, reorganize if-cascades
51 // 20100511 M. Kelsey -- Eliminate three residual random-angles blocks.
52 // 20100511 M. Kelsey -- Bug fix: pi-N two-body final states not correctly
53 // tested for charge-exchange case.
54 // 20100512 M. Kelsey -- Rationalize multiplicity returns to be actual value
55 // 20100512 M. Kelsey -- Add some additional debugging messages for 2-to-2
56 // 20100512 M. Kelsey -- Replace "if (is==)" cascades with switch blocks.
57 // Use G4CascadeInterpolator for angular distributions.
58 // 20100517 M. Kelsey -- Inherit from common base class, make arrays static
59 // 20100519 M. Kelsey -- Use G4InteractionCase to compute "is" values.
60 // 20100625 M. Kelsey -- Two bugs in n-body momentum, last particle recoil
61 // 20100713 M. Kelsey -- Bump collide start message up to verbose > 1
62 // 20100714 M. Kelsey -- Move conservation checking to base class
63 // 20100714 M. Kelsey -- Add sanity check for two-body final state, to ensure
64 // that final state total mass is below etot_scm; also compute
65 // kinematics without "rescaling" (which led to non-conservation)
66 // 20100726 M. Kelsey -- Move remaining std::vector<> buffers to .hh file
67 // 20100804 M. Kelsey -- Add printing of final-state tables, protected by
68 // G4CASCADE_DEBUG_SAMPLER preprocessor flag
69 // 20101019 M. Kelsey -- CoVerity report: check dynamic_cast<> for null
70 // 20110214 M. Kelsey -- Follow G4InuclParticle::Model enumerator migration
71 // 20110718 V. Uzhinskiy -- Drop IL,IM reset for multiplicity-three collisions
72 // 20110718 M. Kelsey -- Use enum names in switch blocks (c.f. G4NucleiModel)
73 // 20110720 M. Kelsey -- Follow interface change for cross-section tables,
74 // eliminating switch blocks.
75 // 20110806 M. Kelsey -- Pre-allocate buffers to avoid memory churn
76 // 20110922 M. Kelsey -- Follow G4InuclParticle::print(ostream&) migration
77 // 20110926 M. Kelsey -- Protect sampleCMcosFor2to2 from unphysical arguments
78 // 20111003 M. Kelsey -- Prepare for gamma-N interactions by checking for
79 // final-state tables instead of particle "isPhoton()"
80 // 20111007 M. Kelsey -- Add gamma-N final-state tables to printFinalState
81 // 20111107 M. Kelsey -- In sampleCMmomentumFor2to2(), hide message about
82 // unrecognized gamma-N initial state behind verbosity.
83 // 20111216 M. Kelsey -- Add diagnostics to generateMomModulesFor2toMany(),
84 // defer allocation of buffer in generateSCMfinalState() so that
85 // multiplicity failures return zero output, and can be trapped.
86 // 20120308 M. Kelsey -- Allow photons to interact with dibaryons (see
87 // changes in G4NucleiModel).
88 // 20120608 M. Kelsey -- Fix variable-name "shadowing" compiler warnings.
89 // 20121206 M. Kelsey -- Add Omega to printFinalStateTables(), remove line
90 // about "preliminary" gamma-N.
91 // 20130129 M. Kelsey -- Add static arrays and interpolators for two-body
92 // angular distributions (addresses MT thread-local issue)
93 // 20130221 M. Kelsey -- Move two-body angular dist classes to factory
94 // 20130306 M. Kelsey -- Use particle-name enums in if-blocks; add comments
95 // to sections of momentum-coefficient matrix; move final state
96 // table printing to G4CascadeChannelTables.
97 // 20130307 M. Kelsey -- Reverse order of dimensions for rmn array
98 // 20130307 M. Kelsey -- Use new momentum generator factory instead of rmn
99 // 20130308 M. Kelsey -- Move 3-body angle calc to G4InuclSpecialFunctions.
100 // 20130422 M. Kelsey -- Move kinematics to G4CascadeFinalStateAlgorithm;
101 // reduce nesting and replicated code in collide().
102 // 20130508 D. Wright -- Implement muon capture
103 // 20130511 M. Kelsey -- Check for neutrinos and skip them in ::collide()
104 // 20141009 M. Kelsey -- Add pion absorption by single nucleons, with
105 // nuclear recoil. Improves pi- capture performance.
106 // 20141030 M. Kelsey -- Check flag for whether to do pi-N absorption
107 // 20141201 M. Kelsey -- Check for null vector return from G4GDecay3
108 // 20141211 M. Kelsey -- Change PIN_ABSORPTION flag to double, probability;
109 // fix handling of boosts for pi-N absorption.
110 // 20150608 M. Kelsey -- Label all while loops as terminating.
111 
113 #include "G4CascadeChannel.hh"
114 #include "G4CascadeChannelTables.hh"
115 #include "G4CascadeParameters.hh"
116 #include "G4CollisionOutput.hh"
117 #include "G4GDecay3.hh"
118 #include "G4InuclParticleNames.hh"
120 #include "G4LorentzConvertor.hh"
122 #include "G4NucleiModel.hh"
123 #include "G4ParticleLargerEkin.hh"
124 #include "G4TwoBodyAngularDist.hh"
125 #include "G4VMultiBodyMomDst.hh"
126 #include "G4VTwoBodyAngDst.hh"
127 #include "Randomize.hh"
128 #include <algorithm>
129 #include <cfloat>
130 #include <vector>
131 
132 using namespace G4InuclParticleNames;
133 using namespace G4InuclSpecialFunctions;
134 
135 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
136 
137 
139  : G4CascadeColliderBase("G4ElementaryParticleCollider"),
140  nucleusA(0), nucleusZ(0) {;}
141 
142 
143 void
146  G4CollisionOutput& output)
147 {
148  if (verboseLevel > 1)
149  G4cout << " >>> G4ElementaryParticleCollider::collide" << G4endl;
150 
151  if (!useEPCollider(bullet,target)) { // Sanity check
152  G4cerr << " ElementaryParticleCollider -> can collide only particle with particle "
153  << G4endl;
154  return;
155  }
156 
157 #ifdef G4CASCADE_DEBUG_SAMPLER
158  static G4bool doPrintTables = true; // Once and only once per job
159  if (doPrintTables) {
161  doPrintTables = false;
162  }
163 #endif
164 
165  interCase.set(bullet, target); // To identify kind of collision
166 
167  if (verboseLevel > 1) G4cout << *bullet << G4endl << *target << G4endl;
168 
169  G4InuclElementaryParticle* particle1 =
170  dynamic_cast<G4InuclElementaryParticle*>(bullet);
171  G4InuclElementaryParticle* particle2 =
172  dynamic_cast<G4InuclElementaryParticle*>(target);
173 
174  if (!particle1 || !particle2) { // Redundant with useEPCollider()
175  G4cerr << " ElementaryParticleCollider -> can only collide hadrons"
176  << G4endl;
177  return;
178  }
179 
180  if (particle1->isNeutrino() || particle2->isNeutrino()) return;
181 
182  // Check for available interaction, or pion+dibaryon special case
184  !particle1->quasi_deutron() && !particle2->quasi_deutron()) {
185  G4cerr << " ElementaryParticleCollider -> cannot collide "
186  << particle1->getDefinition()->GetParticleName() << " with "
187  << particle2->getDefinition()->GetParticleName() << G4endl;
188  return;
189  }
190 
191  G4LorentzConvertor convertToSCM; // Utility to handle frame manipulation
192  if (particle2->nucleon() || particle2->quasi_deutron()) {
193  convertToSCM.setBullet(particle1);
194  convertToSCM.setTarget(particle2);
195  } else {
196  convertToSCM.setBullet(particle2);
197  convertToSCM.setTarget(particle1);
198  }
199 
200  convertToSCM.setVerbose(verboseLevel);
201  convertToSCM.toTheCenterOfMass();
202 
203  G4double etot_scm = convertToSCM.getTotalSCMEnergy();
204 
205  // Generate any particle collision with nucleon
206  if (particle1->nucleon() || particle2->nucleon()) {
207  G4double ekin = convertToSCM.getKinEnergyInTheTRS();
208 
209  // SPECIAL: Very low energy pions may be absorbed by a nucleon
210  if (pionNucleonAbsorption(ekin)) {
211  generateSCMpionNAbsorption(etot_scm, particle1, particle2);
212  } else {
213  generateSCMfinalState(ekin, etot_scm, particle1, particle2);
214  }
215  }
216 
217  // Generate pion or photon collision with quasi-deuteron
218  if (particle1->quasi_deutron() || particle2->quasi_deutron()) {
219  if (!G4NucleiModel::useQuasiDeuteron(particle1->type(),particle2->type()) &&
220  !G4NucleiModel::useQuasiDeuteron(particle2->type(),particle1->type())) {
221  G4cerr << " ElementaryParticleCollider -> can only collide pi,mu,gamma with"
222  << " dibaryons " << G4endl;
223  return;
224  }
225 
226  if (particle1->isMuon() || particle2->isMuon()) {
227  generateSCMmuonAbsorption(etot_scm, particle1, particle2);
228  } else { // Currently, pion absoprtion also handles gammas
229  generateSCMpionAbsorption(etot_scm, particle1, particle2);
230  }
231  }
232 
233  if (particles.empty()) { // No final state possible, pass bullet through
234  if (verboseLevel) {
235  G4cerr << " ElementaryParticleCollider -> failed to collide "
236  << particle1->getMomModule() << " GeV/c "
237  << particle1->getDefinition()->GetParticleName() << " with "
238  << particle2->getDefinition()->GetParticleName() << G4endl;
239  }
240  return;
241  }
242 
243  // Convert final state back to lab frame
244  G4LorentzVector mom; // Buffer to avoid memory churn
245  particleIterator ipart;
246  for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
247  mom = convertToSCM.backToTheLab(ipart->getMomentum());
248  ipart->setMomentum(mom);
249  };
250 
251  // Check conservation in multibody final state
252  if (verboseLevel && !validateOutput(bullet, target, particles)) {
253  G4cout << " incoming particles: \n" << *particle1 << G4endl
254  << *particle2 << G4endl
255  << " outgoing particles: " << G4endl;
256  for(ipart = particles.begin(); ipart != particles.end(); ipart++)
257  G4cout << *ipart << G4endl;
258 
259  G4cout << " <<< Non-conservation in G4ElementaryParticleCollider"
260  << G4endl;
261  }
262 
263  std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());
264  output.addOutgoingParticles(particles);
265 }
266 
267 
268 G4int
269 G4ElementaryParticleCollider::generateMultiplicity(G4int is,
270  G4double ekin) const
271 {
272  G4int mul = 0;
273 
274  const G4CascadeChannel* xsecTable = G4CascadeChannelTables::GetTable(is);
275 
276  if (xsecTable) mul = xsecTable->getMultiplicity(ekin);
277  else {
278  G4cerr << " G4ElementaryParticleCollider: Unknown interaction channel "
279  << is << " - multiplicity not generated " << G4endl;
280  }
281 
282  if(verboseLevel > 3){
283  G4cout << " G4ElementaryParticleCollider::generateMultiplicity: "
284  << " multiplicity = " << mul << G4endl;
285  }
286 
287  return mul;
288 }
289 
290 
291 void
292 G4ElementaryParticleCollider::generateOutgoingPartTypes(G4int is, G4int mult,
293  G4double ekin)
294 {
295  particle_kinds.clear(); // Initialize buffer for generation
296 
297  const G4CascadeChannel* xsecTable = G4CascadeChannelTables::GetTable(is);
298 
299  if (xsecTable)
300  xsecTable->getOutgoingParticleTypes(particle_kinds, mult, ekin);
301  else {
302  G4cerr << " G4ElementaryParticleCollider: Unknown interaction channel "
303  << is << " - outgoing kinds not generated " << G4endl;
304  }
305 
306  return;
307 }
308 
309 
310 void
311 G4ElementaryParticleCollider::generateSCMfinalState(G4double ekin,
312  G4double etot_scm,
313  G4InuclElementaryParticle* particle1,
314  G4InuclElementaryParticle* particle2) {
315  if (verboseLevel > 2) {
316  G4cout << " >>> G4ElementaryParticleCollider::generateSCMfinalState"
317  << G4endl;
318  }
319 
320  fsGenerator.SetVerboseLevel(verboseLevel);
321 
322  const G4int itry_max = 10;
323 
324  G4int type1 = particle1->type();
325  G4int type2 = particle2->type();
326 
327  G4int is = type1 * type2;
328 
329  if (verboseLevel > 3) G4cout << " is " << is << G4endl;
330 
331  G4int multiplicity = 0;
332  G4bool generate = true;
333 
334  G4int itry = 0;
335  while (generate && itry++ < itry_max) { /* Loop checking 08.06.2015 MHK */
336  particles.clear(); // Initialize buffers for this event
337  particle_kinds.clear();
338 
339  // Generate list of final-state particles
340  multiplicity = generateMultiplicity(is, ekin);
341 
342  generateOutgoingPartTypes(is, multiplicity, ekin);
343  if (particle_kinds.empty()) {
344  if (verboseLevel > 3) {
345  G4cout << " generateOutgoingPartTypes failed mult " << multiplicity
346  << G4endl;
347  }
348  continue;
349  }
350 
351  fillOutgoingMasses(); // Fill mass buffer from particle types
352 
353  // Attempt to produce final state kinematics
354  fsGenerator.Configure(particle1, particle2, particle_kinds);
355  generate = !fsGenerator.Generate(etot_scm, masses, scm_momentums);
356  } // while (generate)
357 
358  if (itry >= itry_max) { // Unable to generate valid final state
359  if (verboseLevel > 2)
360  G4cout << " generateSCMfinalState failed " << itry << " attempts"
361  << G4endl;
362  return;
363  }
364 
365  // Store generated momenta into outgoing particles
366 
367  particles.resize(multiplicity); // Preallocate buffer
368  for (G4int i=0; i<multiplicity; i++) {
369  particles[i].fill(scm_momentums[i], particle_kinds[i],
371  }
372 
373  if (verboseLevel > 3) {
374  G4cout << " <<< G4ElementaryParticleCollider::generateSCMfinalState"
375  << G4endl;
376  }
377 
378  return; // Particles buffer filled
379 }
380 
381 
382 // Use generated list of final states to fill mass buffers
383 
384 void G4ElementaryParticleCollider::fillOutgoingMasses() {
385  G4int mult = particle_kinds.size();
386 
387  masses.resize(mult,0.);
388  masses2.resize(mult,0.); // Allows direct [i] setting
389 
390  for (G4int i = 0; i < mult; i++) {
391  masses[i] = G4InuclElementaryParticle::getParticleMass(particle_kinds[i]);
392  masses2[i] = masses[i] * masses[i];
393  }
394 }
395 
396 
397 // generate nucleons momenta for pion or photon absorption by dibaryon
398 // the nucleon distribution assumed to be isotropic in SCM
399 
400 void
401 G4ElementaryParticleCollider::generateSCMpionAbsorption(G4double etot_scm,
402  G4InuclElementaryParticle* particle1,
403  G4InuclElementaryParticle* particle2) {
404  if (verboseLevel > 3)
405  G4cout << " >>> G4ElementaryParticleCollider::generateSCMpionAbsorption"
406  << G4endl;
407 
408  particles.clear(); // Initialize buffers for this event
409  particles.resize(2);
410 
411  particle_kinds.clear();
412 
413  G4int type1 = particle1->type();
414  G4int type2 = particle2->type();
415 
416  // Ensure that absportion is valid (charge conservable)
417  if (!G4NucleiModel::useQuasiDeuteron(type1, type2)) {
418  G4cerr << " pion absorption: "
419  << particle1->getDefinition()->GetParticleName() << " + "
420  << particle2->getDefinition()->GetParticleName() << " -> ?"
421  << G4endl;
422  return;
423  }
424 
425  if (!splitQuasiDeuteron(type2)) return; // Get constituents of [NN]
426 
427  fillOutgoingMasses();
428 
429  G4double a = 0.5 * (etot_scm * etot_scm - masses2[0] - masses2[1]);
430 
431  G4double pmod = std::sqrt((a*a - masses2[0]*masses2[1])
432  / (etot_scm*etot_scm) );
433  G4LorentzVector mom1 = generateWithRandomAngles(pmod, masses[0]);
434  G4LorentzVector mom2;
435  mom2.setVectM(-mom1.vect(), masses[1]);
436 
437  particles[0].fill(mom1, particle_kinds[0], G4InuclParticle::EPCollider);
438  particles[1].fill(mom2, particle_kinds[1], G4InuclParticle::EPCollider);
439 }
440 
441 
442 // generate nucleons momenta for muon absorption by dibaryon
443 // the nucleon distribution assumed to be isotropic in SCM
444 
445 void
446 G4ElementaryParticleCollider::generateSCMmuonAbsorption(G4double etot_scm,
447  G4InuclElementaryParticle* particle1,
448  G4InuclElementaryParticle* particle2)
449 {
450  if (verboseLevel > 3)
451  G4cout << " >>> G4ElementaryParticleCollider::generateSCMmuonAbsorption"
452  << G4endl;
453 
454  // A phase space generator is required for the 3-body final state
455 
456  particles.clear(); // Initialize buffers for this event
457  particles.resize(3);
458 
459  scm_momentums.clear();
460  scm_momentums.resize(3);
461 
462  particle_kinds.clear();
463 
464  G4int type1 = particle1->type();
465  G4int type2 = particle2->type();
466 
467  if (type1 != muonMinus) return; // Sanity check, only mu- absorption
468 
469  // Ensure that absportion is valid (charge conservable)
470  if (!G4NucleiModel::useQuasiDeuteron(type1, type2)) {
471  G4cerr << " mu- absorption: "
472  << particle1->getDefinition()->GetParticleName() << " + "
473  << particle2->getDefinition()->GetParticleName() << " -> ?"
474  << G4endl;
475  return;
476  }
477 
478  if (!splitQuasiDeuteron(type2)) return; // Get constituents of [NN]
479  particle_kinds.push_back(mnu);
480 
481  fillOutgoingMasses();
482 
483  G4GDecay3 breakup(etot_scm, masses[0], masses[1], masses[2]);
484  std::vector<G4ThreeVector> theMomenta = breakup.GetThreeBodyMomenta();
485 
486  if (theMomenta.empty()) {
487  G4cerr << " generateSCMmuonAbsorption: GetThreeBodyMomenta() failed"
488  << " for " << type2 << " dibaryon" << G4endl;
489  particle_kinds.clear();
490  masses.clear();
491  particles.clear();
492  return;
493  }
494 
495  for (size_t i=0; i<3; i++) {
496  scm_momentums[i].setVectM(theMomenta[i], masses[i]);
497  particles[i].fill(scm_momentums[i], particle_kinds[i], G4InuclParticle::EPCollider);
498  }
499 }
500 
501 
502 // generate nucleons momenta for pion absorption by single nucleon
503 
504 void
505 G4ElementaryParticleCollider::generateSCMpionNAbsorption(G4double /*etot_scm*/,
506  G4InuclElementaryParticle* particle1,
507  G4InuclElementaryParticle* particle2) {
508  if (verboseLevel > 3)
509  G4cout << " >>> G4ElementaryParticleCollider::generateSCMpionNAbsorption"
510  << G4endl;
511 
512  particles.clear(); // Initialize buffers for this event
513  particles.resize(1);
514 
515  particle_kinds.clear();
516 
517  G4int type1 = particle1->type();
518  G4int type2 = particle2->type();
519 
520  // Ensure that single-nucleon absportion is valid (charge exchangeable)
521  if ((type1*type2 != pim*pro && type1*type2 != pip*neu)) {
522  G4cerr << " pion-nucleon absorption: "
523  << particle1->getDefinition()->GetParticleName() << " + "
524  << particle2->getDefinition()->GetParticleName() << " -> ?"
525  << G4endl;
526  return;
527  }
528 
529  // Get outgoing nucleon type using charge exchange
530  // Proton code is 1, neutron code is 2, so 3-# swaps them
531  G4int ntype = (particle2->nucleon() ? type2 : type1);
532  G4int outType = 3 - ntype;
533  particle_kinds.push_back(outType);
534 
535  fillOutgoingMasses();
536 
537  // Get mass of residual nucleus (2-ntype = 1 for proton, 0 for neutron)
538  G4double mRecoil =
539  G4InuclNuclei::getNucleiMass(nucleusA-1, nucleusZ-(2-ntype));
540  G4double mRecoil2 = mRecoil*mRecoil;
541 
542  // Recompute Ecm to include nucleus (for recoil kinematics)
543  G4LorentzVector piN4 = particle1->getMomentum() + particle2->getMomentum();
544  G4LorentzVector vsum(0.,0.,0.,mRecoil);
545  vsum += piN4;
546 
547  // Two-body kinematics (nucleon against nucleus) in overall CM system
548  G4double esq_scm = vsum.m2();
549  G4double a = 0.5 * (esq_scm - masses2[0] - mRecoil2);
550 
551  G4double pmod = std::sqrt((a*a - masses2[0]*mRecoil2) / esq_scm );
552  G4LorentzVector mom1 = generateWithRandomAngles(pmod, masses[0]);
553 
554  if (verboseLevel > 3) {
555  G4cout << " outgoing type " << outType << " recoiling on nuclear mass "
556  << mRecoil << "\n a " << a << " p " << pmod << " Ekin "
557  << mom1.e()-masses[0] << G4endl;
558  }
559 
560  mom1.boost(-piN4.boostVector()); // Boost into CM of pi-N collision
561 
562  if (verboseLevel > 3) {
563  G4cout << " in original pi-N frame p(SCM) " << mom1.rho() << " Ekin "
564  << mom1.e()-masses[0] << G4endl;
565  }
566 
567  // Fill only the ejected nucleon
568  particles[0].fill(mom1, particle_kinds[0], G4InuclParticle::EPCollider);
569 }
570 
571 
572 // Evaluate whether interaction is candidate for absorption on nucleon
573 
574 G4bool
575 G4ElementaryParticleCollider::pionNucleonAbsorption(G4double ekin) const {
576  if (verboseLevel > 3)
577  G4cout << " >>> G4ElementaryParticleCollider::pionNucleonAbsorption ?"
578  << " ekin " << ekin << " is " << interCase.hadrons() << G4endl;
579 
580  // Absorption occurs with specified probability
582 
583  // Absorption occurs only for pi- p -> n, or pi+ n -> p
584  // Restrict to "very slow" pions, to allow for some normal scattering
585  return ((interCase.hadrons() == pim*pro ||
586  interCase.hadrons() == pip*neu)
587  && (ekin < 0.05) // 50 MeV kinetic energy or less
588  && (G4UniformRand() < absProb)
589  );
590 }
591 
592 
593 // generate constituents of dibaryon for "explosion"
594 
595 G4bool G4ElementaryParticleCollider::splitQuasiDeuteron(G4int qdtype) {
596  if (qdtype != diproton && qdtype != unboundPN && qdtype != dineutron) {
597  G4cerr << " type " << qdtype << " not dibaryon!" << G4endl;
598  return false;
599  }
600 
601  G4int b2 = qdtype % 10; // Dibaryon codes are 1ab (a=1,2; b=1,2)
602  G4int b1 = (qdtype/10) % 10;
603 
604  particle_kinds.push_back(b1);
605  particle_kinds.push_back(b2);
606 
607  return true;
608 }
Hep3Vector boostVector() const
static const G4CascadeChannel * GetTable(G4int initialState)
G4bool Generate(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
const XML_Char * target
Definition: expat.h:268
G4LorentzVector getMomentum() const
virtual G4bool useEPCollider(G4InuclParticle *bullet, G4InuclParticle *target) const
virtual void getOutgoingParticleTypes(std::vector< G4int > &kinds, G4int mult, G4double ke) const =0
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static G4double piNAbsorption()
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
const G4ParticleDefinition * getDefinition() const
G4int hadrons() const
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
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
void setVectM(const Hep3Vector &spatial, double mass)
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
void setVerbose(G4int vb=0)
bool G4bool
Definition: G4Types.hh:79
HepLorentzVector & boost(double, double, double)
G4double getNucleiMass() const
void Configure(G4InuclElementaryParticle *bullet, G4InuclElementaryParticle *target, const std::vector< G4int > &particle_kinds)
double rho() const
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
G4double getKinEnergyInTheTRS() const
static void Print(std::ostream &os=G4cout)
void addOutgoingParticles(const std::vector< G4InuclElementaryParticle > &particles)
void set(G4InuclParticle *part1, G4InuclParticle *part2)
void generate(const G4double sqrtS, ParticleList &particles)
Generate an event in the CM system.
void SetVerboseLevel(G4int verbose)
double m2() const
virtual G4bool validateOutput(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
std::vector< G4InuclElementaryParticle >::iterator particleIterator
Definition: G4BigBanger.cc:65
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
virtual G4int getMultiplicity(G4double ke) const =0
G4double getMomModule() const
G4GLOB_DLL std::ostream G4cerr
void setTarget(const G4InuclParticle *target)