Geant4_10
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 
106 #include "G4CascadeChannel.hh"
107 #include "G4CascadeChannelTables.hh"
108 #include "G4CollisionOutput.hh"
109 #include "G4GDecay3.hh"
110 #include "G4InuclParticleNames.hh"
112 #include "G4LorentzConvertor.hh"
114 #include "G4NucleiModel.hh"
115 #include "G4ParticleLargerEkin.hh"
116 #include "G4TwoBodyAngularDist.hh"
117 #include "G4VMultiBodyMomDst.hh"
118 #include "G4VTwoBodyAngDst.hh"
119 #include "Randomize.hh"
120 #include <algorithm>
121 #include <cfloat>
122 #include <vector>
123 
124 using namespace G4InuclParticleNames;
125 using namespace G4InuclSpecialFunctions;
126 
127 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
128 
129 
131  : G4CascadeColliderBase("G4ElementaryParticleCollider") {;}
132 
133 
134 void
137  G4CollisionOutput& output)
138 {
139  if (verboseLevel > 1)
140  G4cout << " >>> G4ElementaryParticleCollider::collide" << G4endl;
141 
142  if (!useEPCollider(bullet,target)) { // Sanity check
143  G4cerr << " ElementaryParticleCollider -> can collide only particle with particle "
144  << G4endl;
145  return;
146  }
147 
148 #ifdef G4CASCADE_DEBUG_SAMPLER
149  static G4bool doPrintTables = true; // Once and only once per job
150  if (doPrintTables) {
152  doPrintTables = false;
153  }
154 #endif
155 
156  interCase.set(bullet, target); // To identify kind of collision
157 
158  if (verboseLevel > 1) G4cout << *bullet << G4endl << *target << G4endl;
159 
160  G4InuclElementaryParticle* particle1 =
161  dynamic_cast<G4InuclElementaryParticle*>(bullet);
162  G4InuclElementaryParticle* particle2 =
163  dynamic_cast<G4InuclElementaryParticle*>(target);
164 
165  if (!particle1 || !particle2) { // Redundant with useEPCollider()
166  G4cerr << " ElementaryParticleCollider -> can only collide hadrons"
167  << G4endl;
168  return;
169  }
170 
171  if (particle1->isNeutrino() || particle2->isNeutrino()) return;
172 
173  // Check for available interaction, or pion+dibaryon special case
175  !particle1->quasi_deutron() && !particle2->quasi_deutron()) {
176  G4cerr << " ElementaryParticleCollider -> cannot collide "
177  << particle1->getDefinition()->GetParticleName() << " with "
178  << particle2->getDefinition()->GetParticleName() << G4endl;
179  return;
180  }
181 
182  G4LorentzConvertor convertToSCM; // Utility to handle frame manipulation
183  if (particle2->nucleon() || particle2->quasi_deutron()) {
184  convertToSCM.setBullet(particle1);
185  convertToSCM.setTarget(particle2);
186  } else {
187  convertToSCM.setBullet(particle2);
188  convertToSCM.setTarget(particle1);
189  }
190 
191  convertToSCM.setVerbose(verboseLevel);
192  convertToSCM.toTheCenterOfMass();
193 
194  // Generate any particle collision with nucleon
195  if (particle1->nucleon() || particle2->nucleon()) {
196  G4double ekin = convertToSCM.getKinEnergyInTheTRS();
197  G4double etot_scm = convertToSCM.getTotalSCMEnergy();
198 
199  generateSCMfinalState(ekin, etot_scm, particle1, particle2);
200  }
201 
202  // Generate pion or photon collision with quasi-deuteron
203  if (particle1->quasi_deutron() || particle2->quasi_deutron()) {
204  if (!G4NucleiModel::useQuasiDeuteron(particle1->type(),particle2->type()) &&
205  !G4NucleiModel::useQuasiDeuteron(particle2->type(),particle1->type())) {
206  G4cerr << " ElementaryParticleCollider -> can only collide pi,mu,gamma with"
207  << " dibaryons " << G4endl;
208  return;
209  }
210 
211  G4double etot_scm = convertToSCM.getTotalSCMEnergy();
212 
213  if (particle1->isMuon() || particle2->isMuon()) {
214  generateSCMmuonAbsorption(etot_scm, particle1, particle2);
215  } else { // Currently, pion absoprtion also handles gammas
216  generateSCMpionAbsorption(etot_scm, particle1, particle2);
217  }
218  }
219 
220  if (particles.empty()) { // No final state possible, pass bullet through
221  if (verboseLevel) {
222  G4cerr << " ElementaryParticleCollider -> failed to collide "
223  << particle1->getMomModule() << " GeV/c "
224  << particle1->getDefinition()->GetParticleName() << " with "
225  << particle2->getDefinition()->GetParticleName() << G4endl;
226  }
227  return;
228  }
229 
230  // Convert final state back to lab frame
231  G4LorentzVector mom; // Buffer to avoid memory churn
233  for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
234  mom = convertToSCM.backToTheLab(ipart->getMomentum());
235  ipart->setMomentum(mom);
236  };
237 
238  // Check conservation in multibody final state
239  if (verboseLevel && !validateOutput(bullet, target, particles)) {
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;
245 
246  G4cout << " <<< Non-conservation in G4ElementaryParticleCollider"
247  << G4endl;
248  }
249 
250  std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());
251  output.addOutgoingParticles(particles);
252 }
253 
254 
255 G4int
256 G4ElementaryParticleCollider::generateMultiplicity(G4int is,
257  G4double ekin) const
258 {
259  G4int mul = 0;
260 
261  const G4CascadeChannel* xsecTable = G4CascadeChannelTables::GetTable(is);
262 
263  if (xsecTable) mul = xsecTable->getMultiplicity(ekin);
264  else {
265  G4cerr << " G4ElementaryParticleCollider: Unknown interaction channel "
266  << is << " - multiplicity not generated " << G4endl;
267  }
268 
269  if(verboseLevel > 3){
270  G4cout << " G4ElementaryParticleCollider::generateMultiplicity: "
271  << " multiplicity = " << mul << G4endl;
272  }
273 
274  return mul;
275 }
276 
277 
278 void
279 G4ElementaryParticleCollider::generateOutgoingPartTypes(G4int is, G4int mult,
280  G4double ekin)
281 {
282  particle_kinds.clear(); // Initialize buffer for generation
283 
284  const G4CascadeChannel* xsecTable = G4CascadeChannelTables::GetTable(is);
285 
286  if (xsecTable)
287  xsecTable->getOutgoingParticleTypes(particle_kinds, mult, ekin);
288  else {
289  G4cerr << " G4ElementaryParticleCollider: Unknown interaction channel "
290  << is << " - outgoing kinds not generated " << G4endl;
291  }
292 
293  return;
294 }
295 
296 
297 void
298 G4ElementaryParticleCollider::generateSCMfinalState(G4double ekin,
299  G4double etot_scm,
300  G4InuclElementaryParticle* particle1,
301  G4InuclElementaryParticle* particle2) {
302  if (verboseLevel > 2) {
303  G4cout << " >>> G4ElementaryParticleCollider::generateSCMfinalState"
304  << G4endl;
305  }
306 
307  fsGenerator.SetVerboseLevel(verboseLevel);
308 
309  const G4int itry_max = 10;
310 
311  G4int type1 = particle1->type();
312  G4int type2 = particle2->type();
313 
314  G4int is = type1 * type2;
315 
316  if (verboseLevel > 3) G4cout << " is " << is << G4endl;
317 
318  G4int multiplicity = 0;
319  G4bool generate = true;
320 
321  G4int itry = 0;
322  while (generate && itry++ < itry_max) {
323  particles.clear(); // Initialize buffers for this event
324  particle_kinds.clear();
325 
326  // Generate list of final-state particles
327  multiplicity = generateMultiplicity(is, ekin);
328 
329  generateOutgoingPartTypes(is, multiplicity, ekin);
330  if (particle_kinds.empty()) {
331  if (verboseLevel > 3) {
332  G4cout << " generateOutgoingPartTypes failed mult " << multiplicity
333  << G4endl;
334  }
335  continue;
336  }
337 
338  fillOutgoingMasses(); // Fill mass buffer from particle types
339 
340  // Attempt to produce final state kinematics
341  fsGenerator.Configure(particle1, particle2, particle_kinds);
342  generate = !fsGenerator.Generate(etot_scm, masses, scm_momentums);
343  } // while (generate)
344 
345  if (itry >= itry_max) { // Unable to generate valid final state
346  if (verboseLevel > 2)
347  G4cout << " generateSCMfinalState failed " << itry << " attempts"
348  << G4endl;
349  return;
350  }
351 
352  // Store generated momenta into outgoing particles
353 
354  particles.resize(multiplicity); // Preallocate buffer
355  for (G4int i=0; i<multiplicity; i++) {
356  particles[i].fill(scm_momentums[i], particle_kinds[i],
358  }
359 
360  if (verboseLevel > 3) {
361  G4cout << " <<< G4ElementaryParticleCollider::generateSCMfinalState"
362  << G4endl;
363  }
364 
365  return; // Particles buffer filled
366 }
367 
368 
369 // Use generated list of final states to fill mass buffers
370 
371 void G4ElementaryParticleCollider::fillOutgoingMasses() {
372  G4int mult = particle_kinds.size();
373 
374  masses.resize(mult,0.);
375  masses2.resize(mult,0.); // Allows direct [i] setting
376 
377  for (G4int i = 0; i < mult; i++) {
378  masses[i] = G4InuclElementaryParticle::getParticleMass(particle_kinds[i]);
379  masses2[i] = masses[i] * masses[i];
380  }
381 }
382 
383 
384 // generate nucleons momenta for pion or photon absorption by dibaryon
385 // the nucleon distribution assumed to be isotropic in SCM
386 
387 void
388 G4ElementaryParticleCollider::generateSCMpionAbsorption(G4double etot_scm,
389  G4InuclElementaryParticle* particle1,
390  G4InuclElementaryParticle* particle2) {
391  if (verboseLevel > 3)
392  G4cout << " >>> G4ElementaryParticleCollider::generateSCMpionAbsorption"
393  << G4endl;
394 
395  particles.clear(); // Initialize buffers for this event
396  particles.resize(2);
397 
398  particle_kinds.clear();
399 
400  G4int type1 = particle1->type();
401  G4int type2 = particle2->type();
402 
403  // Ensure that absportion is valid (charge conservable)
404  if (!G4NucleiModel::useQuasiDeuteron(type1, type2)) {
405  G4cerr << " pion absorption: "
406  << particle1->getDefinition()->GetParticleName() << " + "
407  << particle2->getDefinition()->GetParticleName() << " -> ?"
408  << G4endl;
409  return;
410  }
411 
412  if (!splitQuasiDeuteron(type2)) return; // Get constituents of [NN]
413 
414  fillOutgoingMasses();
415 
416  G4double a = 0.5 * (etot_scm * etot_scm - masses2[0] - masses2[1]);
417 
418  G4double pmod = std::sqrt((a*a - masses2[0]*masses2[1])
419  / (masses2[0] + masses2[1] + 2.0*a));
420  G4LorentzVector mom1 = generateWithRandomAngles(pmod, masses[0]);
421  G4LorentzVector mom2;
422  mom2.setVectM(-mom1.vect(), masses[1]);
423 
424  particles[0].fill(mom1, particle_kinds[0], G4InuclParticle::EPCollider);
425  particles[1].fill(mom2, particle_kinds[1], G4InuclParticle::EPCollider);
426 }
427 
428 
429 // generate nucleons momenta for muon absorption by dibaryon
430 // the nucleon distribution assumed to be isotropic in SCM
431 
432 void
433 G4ElementaryParticleCollider::generateSCMmuonAbsorption(G4double etot_scm,
434  G4InuclElementaryParticle* particle1,
435  G4InuclElementaryParticle* particle2)
436 {
437  if (verboseLevel > 3)
438  G4cout << " >>> G4ElementaryParticleCollider::generateSCMmuonAbsorption"
439  << G4endl;
440 
441  // A phase space generator is required for the 3-body final state
442 
443  particles.clear(); // Initialize buffers for this event
444  particles.resize(3);
445 
446  scm_momentums.clear();
447  scm_momentums.resize(3);
448 
449  particle_kinds.clear();
450 
451  G4int type1 = particle1->type();
452  G4int type2 = particle2->type();
453 
454  if (type1 != muonMinus) return; // Sanity check, only mu- absorption
455 
456  // Ensure that absportion is valid (charge conservable)
457  if (!G4NucleiModel::useQuasiDeuteron(type1, type2)) {
458  G4cerr << " mu- absorption: "
459  << particle1->getDefinition()->GetParticleName() << " + "
460  << particle2->getDefinition()->GetParticleName() << " -> ?"
461  << G4endl;
462  return;
463  }
464 
465  if (!splitQuasiDeuteron(type2)) return; // Get constituents of [NN]
466  particle_kinds.push_back(mnu);
467 
468  fillOutgoingMasses();
469 
470  G4GDecay3 breakup(etot_scm, masses[0], masses[1], masses[2]);
471  std::vector<G4ThreeVector> theMomenta = breakup.GetThreeBodyMomenta();
472 
473  for (size_t i=0; i<3; i++) {
474  scm_momentums[i].setVectM(theMomenta[i], masses[i]);
475  particles[i].fill(scm_momentums[i], particle_kinds[i], G4InuclParticle::EPCollider);
476  }
477 }
478 
479 
480 // generate constituents of dibaryon for "explosion"
481 
482 G4bool G4ElementaryParticleCollider::splitQuasiDeuteron(G4int qdtype) {
483  if (qdtype != diproton && qdtype != unboundPN && qdtype != dineutron) {
484  G4cerr << " type " << qdtype << " not dibaryon!" << G4endl;
485  return false;
486  }
487 
488  G4int b2 = qdtype % 10; // Dibaryon codes are 1ab (a=1,2; b=1,2)
489  G4int b1 = (qdtype/10) % 10;
490 
491  particle_kinds.push_back(b1);
492  particle_kinds.push_back(b2);
493 
494  return true;
495 }
static const G4CascadeChannel * GetTable(G4int initialState)
G4bool Generate(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
tuple a
Definition: test.py:11
virtual G4bool useEPCollider(G4InuclParticle *bullet, G4InuclParticle *target) const
virtual void getOutgoingParticleTypes(std::vector< G4int > &kinds, G4int mult, G4double ke) const =0
Int_t ipart
Definition: Style.C:10
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
G4int hadrons() const
subroutine sort(A, N)
Definition: dpm25nuc7.f:4670
static G4bool useQuasiDeuteron(G4int ptype, G4int qdtype=0)
const XML_Char * target
Definition: expat.h:268
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
G4GLOB_DLL std::ostream G4cout
void setVerbose(G4int vb=0)
bool G4bool
Definition: G4Types.hh:79
void Configure(G4InuclElementaryParticle *bullet, G4InuclElementaryParticle *target, const std::vector< G4int > &particle_kinds)
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 SetVerboseLevel(G4int verbose)
virtual G4bool validateOutput(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
std::vector< G4InuclElementaryParticle >::iterator particleIterator
Definition: G4BigBanger.cc:64
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
virtual G4int getMultiplicity(G4double ke) const =0
G4ParticleDefinition * getDefinition() const
G4double getMomModule() const
G4GLOB_DLL std::ostream G4cerr
void setTarget(const G4InuclParticle *target)