Context Navigation

source: trunk/source/processes/hadronic/models/cascade/cascade/src/G4ElementaryParticleCollider.cc @ 1315

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago
update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 38.3 KB

Line
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	// $Id: G4ElementaryParticleCollider.cc,v 1.64 2010/05/21 18:26:16 mkelsey Exp $
26	// Geant4 tag: $Name: geant4-09-04-beta-cand-01 $
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
61	#include "G4ElementaryParticleCollider.hh"
62
63	#include "G4CascadePiMinusNChannel.hh"
64	#include "G4CascadePiMinusPChannel.hh"
65	#include "G4CascadePiPlusNChannel.hh"
66	#include "G4CascadePiPlusPChannel.hh"
67	#include "G4CascadePiZeroNChannel.hh"
68	#include "G4CascadePiZeroPChannel.hh"
69	#include "G4CascadeKplusPChannel.hh"
70	#include "G4CascadeKplusNChannel.hh"
71	#include "G4CascadeKzeroPChannel.hh"
72	#include "G4CascadeKzeroNChannel.hh"
73	#include "G4CascadeKminusPChannel.hh"
74	#include "G4CascadeKminusNChannel.hh"
75	#include "G4CascadeKzeroBarPChannel.hh"
76	#include "G4CascadeKzeroBarNChannel.hh"
77	#include "G4CascadeNNChannel.hh"
78	#include "G4CascadeNPChannel.hh"
79	#include "G4CascadePPChannel.hh"
80	#include "G4CascadeLambdaPChannel.hh"
81	#include "G4CascadeLambdaNChannel.hh"
82	#include "G4CascadeSigmaPlusPChannel.hh"
83	#include "G4CascadeSigmaPlusNChannel.hh"
84	#include "G4CascadeSigmaZeroPChannel.hh"
85	#include "G4CascadeSigmaZeroNChannel.hh"
86	#include "G4CascadeSigmaMinusPChannel.hh"
87	#include "G4CascadeSigmaMinusNChannel.hh"
88	#include "G4CascadeXiZeroPChannel.hh"
89	#include "G4CascadeXiZeroNChannel.hh"
90	#include "G4CascadeXiMinusPChannel.hh"
91	#include "G4CascadeXiMinusNChannel.hh"
92
93	#include "G4CascadeInterpolator.hh"
94	#include "G4CollisionOutput.hh"
95	#include "G4InuclParticleNames.hh"
96	#include "G4InuclSpecialFunctions.hh"
97	#include "G4ParticleLargerEkin.hh"
98	#include "G4LorentzConvertor.hh"
99	#include "Randomize.hh"
100	#include <algorithm>
101	#include <vector>
102
103	using namespace G4InuclParticleNames;
104	using namespace G4InuclSpecialFunctions;
105
106	typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
107
108
109	G4ElementaryParticleCollider::G4ElementaryParticleCollider()
110	: G4VCascadeCollider("G4ElementaryParticleCollider") {}
111
112
113	void
114	G4ElementaryParticleCollider::collide(G4InuclParticle* bullet,
115	G4InuclParticle* target,
116	G4CollisionOutput& output)
117	{
118	if (!useEPCollider(bullet,target)) { // Sanity check
119	G4cerr << " ElementaryParticleCollider -> can collide only particle with particle "
120	<< G4endl;
121	return;
122	}
123
124	interCase.set(bullet, target); // To identify kind of collision
125
126	G4InuclElementaryParticle* particle1 =
127	dynamic_cast<G4InuclElementaryParticle*>(bullet);
128	G4InuclElementaryParticle* particle2 =
129	dynamic_cast<G4InuclElementaryParticle*>(target);
130
131	if (particle1->isPhoton() \|\| particle2->isPhoton()) {
132	G4cout << " ElementaryParticleCollider -> cannot collide photons "
133	<< G4endl;
134	} else {
135
136	// Generate nucleon or pion collision with nucleon
137	// or pion with quasi-deuteron
138
139	if (particle1->nucleon() \|\| particle2->nucleon()) { // ok
140	G4LorentzConvertor convertToSCM;
141	if(particle2->nucleon()) {
142	convertToSCM.setBullet(particle1);
143	convertToSCM.setTarget(particle2);
144	} else {
145	convertToSCM.setBullet(particle2);
146	convertToSCM.setTarget(particle1);
147	};
148	convertToSCM.toTheCenterOfMass();
149	G4double ekin = convertToSCM.getKinEnergyInTheTRS();
150	G4double etot_scm = convertToSCM.getTotalSCMEnergy();
151	G4double pscm = convertToSCM.getSCMMomentum();
152
153	generateSCMfinalState(ekin, etot_scm, pscm, particle1, particle2,
154	&convertToSCM);
155
156	if(verboseLevel > 2){
157	G4cout << " particles " << particles.size() << G4endl;
158
159	for(G4int i = 0; i < G4int(particles.size()); i++)
160	particles[i].printParticle();
161
162	}
163	if(!particles.empty()) { // convert back to Lab
164	G4LorentzVector mom; // Buffer to avoid memory churn
165	particleIterator ipart;
166	for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
167	mom = convertToSCM.backToTheLab(ipart->getMomentum());
168	ipart->setMomentum(mom);
169	};
170	std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());
171	output.addOutgoingParticles(particles);
172	};
173	} else {
174	if(particle1->quasi_deutron() \|\| particle2->quasi_deutron()) {
175	if(particle1->pion() \|\| particle2->pion()) {
176	G4LorentzConvertor convertToSCM;
177	if(particle1->pion()) {
178	convertToSCM.setBullet(particle1);
179	convertToSCM.setTarget(particle2);
180	} else {
181	convertToSCM.setBullet(particle2);
182	convertToSCM.setTarget(particle1);
183	};
184	convertToSCM.toTheCenterOfMass();
185	G4double etot_scm = convertToSCM.getTotalSCMEnergy();
186
187	generateSCMpionAbsorption(etot_scm, particle1, particle2);
188
189	if(!particles.empty()) { // convert back to Lab
190	G4LorentzVector mom; // Buffer to avoid memory churn
191	particleIterator ipart;
192	for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
193	mom = convertToSCM.backToTheLab(ipart->getMomentum());
194	ipart->setMomentum(mom);
195	};
196	std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());
197	output.addOutgoingParticles(particles);
198	};
199
200	} else {
201	G4cout << " ElementaryParticleCollider -> can only collide pions with deuterons "
202	<< G4endl;
203	};
204	} else {
205	G4cout << " ElementaryParticleCollider -> can only collide something with nucleon or deuteron "
206	<< G4endl;
207	};
208	};
209
210	}
211
212	}
213
214
215	G4int
216	G4ElementaryParticleCollider::generateMultiplicity(G4int is,
217	G4double ekin) const
218	{
219	G4int mul = 0;
220
221	switch (is) {
222	case 1: mul = G4CascadePPChannel::getMultiplicity(ekin); break;
223	case 2: mul = G4CascadeNPChannel::getMultiplicity(ekin); break;
224	case 3: mul = G4CascadePiPlusPChannel::getMultiplicity(ekin); break;
225	case 4: mul = G4CascadeNNChannel::getMultiplicity(ekin); break;
226	case 5: mul = G4CascadePiMinusPChannel::getMultiplicity(ekin); break;
227	case 6: mul = G4CascadePiPlusNChannel::getMultiplicity(ekin); break;
228	case 7: mul = G4CascadePiZeroPChannel::getMultiplicity(ekin); break;
229	case 10: mul = G4CascadePiMinusNChannel::getMultiplicity(ekin); break;
230	case 11: mul = G4CascadeKplusPChannel::getMultiplicity(ekin); break;
231	case 13: mul = G4CascadeKminusPChannel::getMultiplicity(ekin); break;
232	case 14: mul = G4CascadePiZeroNChannel::getMultiplicity(ekin); break;
233	case 15: mul = G4CascadeKzeroPChannel::getMultiplicity(ekin); break;
234	case 17: mul = G4CascadeKzeroBarPChannel::getMultiplicity(ekin); break;
235	case 21: mul = G4CascadeLambdaPChannel::getMultiplicity(ekin); break;
236	case 22: mul = G4CascadeKplusNChannel::getMultiplicity(ekin); break;
237	case 23: mul = G4CascadeSigmaPlusPChannel::getMultiplicity(ekin); break;
238	case 25: mul = G4CascadeSigmaZeroPChannel::getMultiplicity(ekin); break;
239	case 26: mul = G4CascadeKminusNChannel::getMultiplicity(ekin); break;
240	case 27: mul = G4CascadeSigmaMinusPChannel::getMultiplicity(ekin); break;
241	case 29: mul = G4CascadeXiZeroPChannel::getMultiplicity(ekin); break;
242	case 30: mul = G4CascadeKzeroNChannel::getMultiplicity(ekin); break;
243	case 31: mul = G4CascadeXiMinusPChannel::getMultiplicity(ekin); break;
244	case 34: mul = G4CascadeKzeroBarNChannel::getMultiplicity(ekin); break;
245	case 42: mul = G4CascadeLambdaNChannel::getMultiplicity(ekin); break;
246	case 46: mul = G4CascadeSigmaPlusNChannel::getMultiplicity(ekin); break;
247	case 50: mul = G4CascadeSigmaZeroNChannel::getMultiplicity(ekin); break;
248	case 54: mul = G4CascadeSigmaMinusNChannel::getMultiplicity(ekin); break;
249	case 58: mul = G4CascadeXiZeroNChannel::getMultiplicity(ekin); break;
250	case 62: mul = G4CascadeXiMinusNChannel::getMultiplicity(ekin); break;
251	default:
252	G4cerr << " G4ElementaryParticleCollider: Unknown interaction channel "
253	<< is << " - multiplicity not generated " << G4endl;
254	}
255
256	if(verboseLevel > 3){
257	G4cout << " G4ElementaryParticleCollider::generateMultiplicity: "
258	<< " multiplicity = " << mul << G4endl;
259	}
260
261	return mul;
262	}
263
264
265	void
266	G4ElementaryParticleCollider::generateSCMfinalState(G4double ekin,
267	G4double etot_scm,
268	G4double pscm,
269	G4InuclElementaryParticle* particle1,
270	G4InuclElementaryParticle* particle2,
271	G4LorentzConvertor* toSCM) {
272	if (verboseLevel > 3) {
273	G4cout << " >>> G4ElementaryParticleCollider::generateSCMfinalState"
274	<< G4endl;
275	}
276
277	const G4double ang_cut = 0.9999;
278	const G4double difr_const = 0.3678794;
279	const G4int itry_max = 10;
280	G4InuclElementaryParticle dummy;
281
282	G4int type1 = particle1->type();
283	G4int type2 = particle2->type();
284
285	G4int is = type1 * type2;
286
287	if(verboseLevel > 3){
288	G4cout << " is " << is << G4endl;
289	}
290
291	G4int multiplicity = 0;
292	G4bool generate = true;
293
294	particles.clear(); // Initialize buffers for this event
295	particle_kinds.clear();
296
297	while (generate) {
298	if(multiplicity == 0) {
299	multiplicity = generateMultiplicity(is, ekin);
300	} else {
301	multiplicity = generateMultiplicity(is, ekin);
302	particle_kinds.clear();
303	}
304
305	// Generate list of final-state particles
306	generateOutgoingPartTypes(is, multiplicity, ekin);
307
308	if (particle_kinds.empty()) continue;
309	// G4cout << " Particle kinds = " ;
310	// for (G4int i = 0; i < multiplicity; i++) G4cout << particle_kinds[i] << " , " ;
311	// G4cout << G4endl;
312
313	if (multiplicity == 2) {
314	// Identify charge or strangeness exchange (non-elastic scatter)
315	G4int finaltype = particle_kinds[0]*particle_kinds[1];
316	G4int kw = (finaltype != is) ? 2 : 1;
317
318	G4double pmod = pscm; // May need to rescale momentum
319	if (kw == 2) {
320	G4double m1 = dummy.getParticleMass(particle_kinds[0]);
321	m1 *= m1;
322	G4double m2 = dummy.getParticleMass(particle_kinds[1]);
323	m2 *= m2;
324	G4double a = 0.5 * (etot_scm * etot_scm - m1 - m2);
325	G4double em = a * a - m1 * m2;
326
327	if (em > 0) { // It is possible to rescale
328	pmod = std::sqrt( em / (m1 + m2 + 2.0 * a));
329	}
330	}
331
332	G4LorentzVector mom = sampleCMmomentumFor2to2(is, kw, ekin, pmod);
333
334	if (verboseLevel > 3){
335	G4cout << " Particle kinds = " << particle_kinds[0] << " , "
336	<< particle_kinds[1] << G4endl
337	<< " before rotation px " << mom.x() << " py " << mom.y()
338	<< " pz " << mom.z() << G4endl;
339	}
340
341	mom = toSCM->rotate(mom);
342
343	if (verboseLevel > 3){
344	G4cout << " after rotation px " << mom.x() << " py " << mom.y() <<
345	" pz " << mom.z() << G4endl;
346	}
347	G4LorentzVector mom1 = -mom;
348
349	particles.push_back(G4InuclElementaryParticle(mom, particle_kinds[0], 3));
350	particles.push_back(G4InuclElementaryParticle(mom1, particle_kinds[1],3));
351	generate = false;
352
353	} else { // 2 -> many
354	G4int itry = 0;
355	G4bool bad = true;
356	G4int knd_last = particle_kinds[multiplicity - 1];
357	G4double mass_last = dummy.getParticleMass(knd_last);
358
359	if (verboseLevel > 3){
360	G4cout << " knd_last " << knd_last << " mass " << mass_last << G4endl;
361	}
362
363	while (bad && itry < itry_max) {
364	itry++;
365
366	if (verboseLevel > 3){
367	G4cout << " itry in while " << itry << G4endl;
368	}
369
370	generateMomModules(multiplicity, is, ekin, etot_scm);
371	if (G4int(modules.size()) != multiplicity) continue;
372
373	if (multiplicity == 3) {
374	G4LorentzVector mom3 =
375	particleSCMmomentumFor2to3(is, knd_last, ekin, modules[2]);
376
377	mom3 = toSCM->rotate(mom3);
378
379	// generate the momentum of first
380	G4double ct = -0.5 * (modules[2] * modules[2] +
381	modules[0] * modules[0] -
382	modules[1] * modules[1]) /
383	modules[2] / modules[0];
384
385	if(std::fabs(ct) < ang_cut) {
386
387	if(verboseLevel > 2){
388	G4cout << " ok for mult " << multiplicity << G4endl;
389	}
390
391	G4LorentzVector mom1 = generateWithFixedTheta(ct, modules[0]);
392
393	mom1 = toSCM->rotate(mom3, mom1);
394
395	G4LorentzVector mom2 = -(mom3 + mom1);
396
397	bad = false;
398	generate = false;
399
400	particles.push_back(G4InuclElementaryParticle(mom1, particle_kinds[0],3));
401	particles.push_back(G4InuclElementaryParticle(mom2, particle_kinds[1],3));
402	particles.push_back(G4InuclElementaryParticle(mom3, particle_kinds[2],3));
403	};
404	} else { // multiplicity > 3
405	// generate first mult - 2 momentums
406	std::vector<G4LorentzVector> scm_momentums;
407	G4LorentzVector tot_mom;
408
409	for (G4int i = 0; i < multiplicity - 2; i++) {
410	G4double p0 = particle_kinds[i] < 3 ? 0.36 : 0.25;
411	G4double alf = 1.0 / p0 / (p0 - (modules[i] + p0) *
412	std::exp(-modules[i] / p0));
413	G4double st = 2.0;
414	G4int itry1 = 0;
415
416	while (std::fabs(st) > ang_cut && itry1 < itry_max) {
417	itry1++;
418	G4double s1 = modules[i] * inuclRndm();
419	G4double s2 = alf * difr_const * p0 * inuclRndm();
420
421	if(verboseLevel > 3){
422	G4cout << " s1 * alf * std::exp(-s1 / p0) "
423	<< s1 * alf * std::exp(-s1 / p0)
424	<< " s2 " << s2 << G4endl;
425	}
426
427	if(s1 * alf * std::exp(-s1 / p0) > s2) st = s1 / modules[i];
428
429	};
430
431	if(verboseLevel > 3){
432	G4cout << " itry1 " << itry1 << " i " << i << " st " << st
433	<< G4endl;
434	}
435
436	if(itry1 == itry_max) {
437	if(verboseLevel > 2){
438	G4cout << " high energy angles generation: itry1 " << itry1
439	<< G4endl;
440	}
441
442	st = 0.5 * inuclRndm();
443	};
444
445	G4double ct = std::sqrt(1.0 - st * st);
446	if (inuclRndm() > 0.5) ct = -ct;
447
448	G4LorentzVector mom = generateWithFixedTheta(ct,modules[i]);
449
450	tot_mom += mom;
451
452	scm_momentums.push_back(mom);
453	};
454
455	// handle last two
456	G4double tot_mod = tot_mom.rho();
457	G4double ct = -0.5 * (tot_mod * tot_mod +
458	modules[multiplicity - 2] * modules[multiplicity - 2] -
459	modules[multiplicity - 1] * modules[multiplicity - 1]) / tot_mod /
460	modules[multiplicity - 2];
461
462	if (verboseLevel > 2){
463	G4cout << " ct last " << ct << G4endl;
464	}
465
466	if (std::fabs(ct) < ang_cut) {
467
468	G4int i(0);
469	for (i = 0; i < multiplicity - 2; i++)
470	scm_momentums[i] = toSCM->rotate(scm_momentums[i]);
471
472	tot_mom = toSCM->rotate(tot_mom);
473
474	G4LorentzVector mom =
475	generateWithFixedTheta(ct, modules[multiplicity - 2]);
476
477	mom = toSCM->rotate(tot_mom, mom);
478	scm_momentums.push_back(mom);
479
480	// and the last one
481	G4LorentzVector mom1 = -(mom+tot_mom);
482
483	scm_momentums.push_back(mom1);
484	bad = false;
485	generate = false;
486
487	if (verboseLevel > 2){
488	G4cout << " ok for mult " << multiplicity << G4endl;
489	}
490
491	for (i = 0; i < multiplicity; i++) {
492	particles.push_back(G4InuclElementaryParticle(scm_momentums[i], particle_kinds[i],3));
493	};
494	};
495	};
496	};
497
498	if (itry == itry_max) {
499	if (verboseLevel > 2){
500	G4cout << " cannot generate the distr. for mult " << multiplicity <<
501	G4endl << " and set it to " << multiplicity - 1 << G4endl;
502	}
503	};
504	};
505	};
506
507	if (verboseLevel > 3) {
508	G4cout << " <<< G4ElementaryParticleCollider::generateSCMfinalState"
509	<< G4endl;
510	}
511
512	return; // Particles buffer filled
513	}
514
515	void
516	G4ElementaryParticleCollider::generateMomModules(G4int mult,
517	G4int is,
518	G4double ekin,
519	G4double etot_cm) {
520	if (verboseLevel > 3) {
521	G4cout << " >>> G4ElementaryParticleCollider::generateMomModules"
522	<< G4endl;
523	}
524
525	if (verboseLevel > 2){
526	G4cout << " mult " << mult << " is " << is << " ekin " << ekin
527	<< " etot_cm " << etot_cm << G4endl;
528	}
529
530	const G4int itry_max = 10;
531	const G4double small = 1.e-10;
532	G4InuclElementaryParticle dummy;
533	G4int itry = 0;
534
535	// FIXME: Code below wants to set modules[i] directly. Bad practice
536	modules.clear(); // Initialize buffer for this attempt
537	modules.insert(modules.begin(), mult, 0.);
538
539	std::vector<G4double> masses2(mult);
540
541	for (G4int i = 0; i < mult; i++) {
542	G4double mass = dummy.getParticleMass(particle_kinds[i]);
543	masses2[i] = mass * mass;
544	};
545
546	G4double mass_last = std::sqrt(masses2[mult - 1]);
547
548	if (verboseLevel > 3){
549	G4cout << " knd_last " << particle_kinds[mult - 1] << " mlast "
550	<< mass_last << G4endl;
551	}
552
553	while (itry < itry_max) {
554	itry++;
555	if(verboseLevel > 3){
556	G4cout << " itry in generateMomModules " << itry << G4endl;
557	}
558
559	G4int ilast = -1;
560	G4double eleft = etot_cm;
561
562	for (G4int i = 0; i < mult - 1; i++) {
563
564	G4double pmod =
565	getMomModuleFor2toMany(is, mult, particle_kinds[i], ekin);
566
567	if (pmod < small) break;
568	eleft -= std::sqrt(pmod * pmod + masses2[i]);
569
570	if (verboseLevel > 3){
571	G4cout << " kp " << particle_kinds[i] << " pmod " << pmod << " mass2 "
572	<< masses2[i] << G4endl;
573	G4cout << " x1 " << eleft - mass_last << G4endl;
574	}
575
576	if (eleft <= mass_last) break;
577	ilast++;
578	modules[i] = pmod;
579	};
580
581	if (ilast == mult - 2) {
582	G4double plast = eleft * eleft - masses2[mult - 1];
583	if (verboseLevel > 2){
584	G4cout << " plast ** 2 " << plast << G4endl;
585	}
586
587	if (plast > small) {
588	plast = std::sqrt(plast);
589	modules[mult - 1] = plast;
590
591	if (mult == 3) {
592	if (satisfyTriangle(modules)) return;
593	} else return;
594	}
595	}
596	}
597
598	modules.clear(); // Something went wrong, throw away partial
599	return;
600	}
601
602
603	G4bool
604	G4ElementaryParticleCollider::satisfyTriangle(
605	const std::vector<G4double>& modules) const
606	{
607	if (verboseLevel > 3) {
608	G4cout << " >>> G4ElementaryParticleCollider::satisfyTriangle"
609	<< G4endl;
610	}
611
612	G4bool good = ( (modules.size() != 3) \|\|
613	!(std::fabs(modules[1] - modules[2]) > modules[0] \|\|
614	modules[0] > modules[1] + modules[2] \|\|
615	std::fabs(modules[0] - modules[2]) > modules[1] \|\|
616	modules[1] > modules[0] + modules[2] \|\|
617	std::fabs(modules[0] - modules[1]) > modules[2] \|\|
618	modules[2] > modules[1] + modules[0]));
619
620	return good;
621	}
622
623
624	void
625	G4ElementaryParticleCollider::generateOutgoingPartTypes(G4int is, G4int mult,
626	G4double ekin)
627	{
628	particle_kinds.clear(); // Initialize buffer for generation
629
630	switch (is) {
631	case 1: G4CascadePPChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
632	case 2: G4CascadeNPChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
633	case 3: G4CascadePiPlusPChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
634	case 4: G4CascadeNNChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
635	case 5: G4CascadePiMinusPChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
636	case 6: G4CascadePiPlusNChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
637	case 7: G4CascadePiZeroPChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
638	case 10: G4CascadePiMinusNChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
639	case 11: G4CascadeKplusPChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
640	case 13: G4CascadeKminusPChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
641	case 14: G4CascadePiZeroNChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
642	case 15: G4CascadeKzeroPChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
643	case 17: G4CascadeKzeroBarPChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
644	case 21: G4CascadeLambdaPChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
645	case 22: G4CascadeKplusNChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
646	case 23: G4CascadeSigmaPlusPChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
647	case 25: G4CascadeSigmaZeroPChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
648	case 26: G4CascadeKminusNChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
649	case 27: G4CascadeSigmaMinusPChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
650	case 29: G4CascadeXiZeroPChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
651	case 30: G4CascadeKzeroNChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
652	case 31: G4CascadeXiMinusPChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
653	case 34: G4CascadeKzeroBarNChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
654	case 42: G4CascadeLambdaNChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
655	case 46: G4CascadeSigmaPlusNChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
656	case 50: G4CascadeSigmaZeroNChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
657	case 54: G4CascadeSigmaMinusNChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
658	case 58: G4CascadeXiZeroNChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
659	case 62: G4CascadeXiMinusNChannel::getOutgoingParticleTypes(particle_kinds, mult, ekin); break;
660	default:
661	G4cout << " G4ElementaryParticleCollider: Unknown interaction channel "
662	<< is << " - outgoing kinds not generated " << G4endl;
663	}
664
665	return;
666	}
667
668
669	G4double
670	G4ElementaryParticleCollider::getMomModuleFor2toMany(G4int is, G4int mult,
671	G4int knd,
672	G4double ekin) const
673	{
674	if (verboseLevel > 3) {
675	G4cout << " >>> G4ElementaryParticleCollider::getMomModuleFor2toMany"
676	<< G4endl;
677	}
678
679	G4double S = inuclRndm();
680	G4double PS = 0.0;
681	G4double PR = 0.0;
682	G4double PQ = 0.0;
683	G4int KM = 2;
684	G4int IL = 4;
685	G4int JK = 4;
686	G4int JM = 2;
687	G4int IM = 3;
688
689	if(is == 1 \|\| is == 2 \|\| is == 4) KM = 1;
690	if(mult == 3) { IM = 0; IL = 0; };
691	if(knd == 1 \|\| knd == 2) JK = 0;
692
693	for(G4int i = 0; i < 4; i++) {
694	G4double V = 0.0;
695	for(G4int k = 0; k < 4; k++) V += rmn[k + JK][i + IL][KM - 1] * std::pow(ekin, k);
696	PR += V * std::pow(S, i);
697	PQ += V;
698	}
699
700	if(knd == 1 \|\| knd == 2) JM = 1;
701	for(G4int m = 0; m < 3; m++) PS += rmn[8 + IM + m][7 + JM][KM - 1] * std::pow(ekin, m);
702	G4double PRA = PS * std::sqrt(S) * (PR + (1 - PQ) * (SSS*S));
703
704	return std::fabs(PRA);
705	}
706
707
708	G4LorentzVector
709	G4ElementaryParticleCollider::particleSCMmomentumFor2to3(
710	G4int is,
711	G4int knd,
712	G4double ekin,
713	G4double pmod) const
714	{
715	if (verboseLevel > 3) {
716	G4cout << " >>> G4ElementaryParticleCollider::particleSCMmomentumFor2to3"
717	<< G4endl;
718	}
719
720	const G4int itry_max = 100;
721	G4double ct = 2.0;
722	G4int K = 3;
723	G4int J = 1;
724
725	if(is == 1 \|\| is == 2 \|\| is == 4) K = 1;
726
727	if(knd == 1 \|\| knd == 2) J = 0;
728
729	G4int itry = 0;
730
731	while(std::fabs(ct) > 1.0 && itry < itry_max) {
732	itry++;
733	G4double S = inuclRndm();
734	G4double U = 0.0;
735	G4double W = 0.0;
736
737	for(G4int l = 0; l < 4; l++) {
738	G4double V = 0.0;
739
740	for(G4int m = 0; m < 4; m++) {
741	V += abn[m][l][K + J - 1] * std::pow(ekin, m);
742	};
743
744	U += V;
745	W += V * std::pow(S, l);
746	};
747	ct = 2.0 * std::sqrt(S) * (W + (1.0 - U) * (SSS*S)) - 1.0;
748	};
749
750	if(itry == itry_max) {
751
752	if(verboseLevel > 2){
753	G4cout << " particleSCMmomentumFor2to3 -> itry = itry_max " << itry << G4endl;
754	}
755
756	ct = 2.0 * inuclRndm() - 1.0;
757	};
758
759	return generateWithFixedTheta(ct, pmod);
760	}
761
762
763	G4LorentzVector
764	G4ElementaryParticleCollider::sampleCMmomentumFor2to2(G4int is, G4int kw,
765	G4double ekin,
766	G4double pscm) const
767	{
768	if (verboseLevel > 3)
769	G4cout << " >>> G4ElementaryParticleCollider::sampleCMmomentumFor2to2"
770	<< " is " << is << " kw " << kw << " ekin " << ekin << " pscm "
771	<< pscm << G4endl;
772
773	G4double pA=0.0, pC=0.0, pCos=0.0, pFrac=0.0; // Angular parameters
774
775	// Arrays below are parameters for two-exponential sampling of angular
776	// distributions of two-body scattering in the specified channels
777
778	if (is == 1 \|\| is == 2 \|\| is == 4 \|\|
779	is == 21 \|\| is == 23 \|\| is == 25 \|\| is == 27 \|\| is ==29 \|\| is == 31 \|\|
780	is == 42 \|\| is == 46 \|\| is == 50 \|\| is == 54 \|\| is ==58 \|\| is == 62) {
781	// nucleon-nucleon or hyperon-nucleon
782	if (verboseLevel > 3) G4cout << " nucleon/hyperon elastic" << G4endl;
783
784	static const G4double nnke[9] = { 0.0, 0.44, 0.59, 0.80, 1.00, 2.24, 4.40, 6.15, 10.00};
785	static const G4double nnA[9] = { 0.0, 0.34, 2.51, 4.59, 4.2, 5.61, 6.38, 7.93, 8.7};
786	static const G4double nnC[9] = { 0.0, 0.0, 1.21, 1.54, 1.88, 1.24, 1.91, 4.04, 8.7};
787	static const G4double nnCos[9] = {-1.0, -1.0, 0.4226, 0.4226, 0.4384, 0.7193, 0.8788, 0.9164, 0.95};
788	static const G4double nnFrac[9] = {1.0, 1.0, 0.4898, 0.7243, 0.7990, 0.8892, 0.8493, 0.9583, 1.0};
789
790	static G4CascadeInterpolator<9> interp(nnke); // Only need one!
791	pA = interp.interpolate(ekin, nnA);
792	pC = interp.interpolate(ekin, nnC);
793	pCos = interp.interpolate(ekin, nnCos);
794	pFrac = interp.interpolate(ekin, nnFrac);
795	} else if (kw == 2) {
796	// pion charge exchange (pi-p -> pi0n, pi+n -> pi0p, pi0p -> pi+n, pi0n -> pi-p
797	if (verboseLevel > 3) G4cout << " pion-nucleon charge exchange " << G4endl;
798
799	static const G4double nnke[10] = {0.0, 0.062, 0.12, 0.217, 0.533, 0.873, 1.34, 2.86, 5.86, 10.0};
800	static const G4double nnA[10] = {0.0, 0.0, 2.48, 7.93, 10.0, 9.78, 5.08, 8.13, 8.13, 8.13};
801	static const G4double nnC[10] = {0.0, -39.58, -12.55, -4.38, 1.81, -1.99, -0.33, 1.2, 1.43, 8.13};
802	static const G4double nnCos[10] = {1.0, 1.0, 0.604, -0.033, 0.25, 0.55, 0.65, 0.80, 0.916, 0.916};
803	static const G4double nnFrac[10] = {0.0, 0.0, 0.1156, 0.5832, 0.8125, 0.3357, 0.3269, 0.7765, 0.8633, 1.0};
804
805	static G4CascadeInterpolator<10> interp(nnke); // Only need one!
806	pA = interp.interpolate(ekin, nnA);
807	pC = interp.interpolate(ekin, nnC);
808	pCos = interp.interpolate(ekin, nnCos);
809	pFrac = interp.interpolate(ekin, nnFrac);
810	} else if (is == 3 \|\| is == 7 \|\| is == 14 \|\| is == 10 \|\| is == 11 \|\|
811	is == 30 \|\| is == 17 \|\| is == 26) {
812	// pi+p, pi0p, pi0n, pi-n, k+p, k0n, k0bp, or k-n
813	if (verboseLevel > 3) G4cout << " meson-nucleon elastic (1)" << G4endl;
814
815	static const G4double nnke[10] = {0.0, 0.062, 0.12, 0.217, 0.533, 0.873, 1.34, 2.86, 5.86, 10.0};
816	static const G4double nnA[10] = {0.0, 0.0, 27.58, 19.83, 6.46, 4.59, 6.47, 6.68, 6.43, 6.7};
817	static const G4double nnC[10] = {0.0, -26.4, -30.55, -19.42, -5.05, -5.24, -1.00, 2.14, 2.9, 6.7};
818	static const G4double nnCos[10] = {1.0, 1.0, 0.174, -0.174, -0.7, -0.295, 0.5, 0.732, 0.837, 0.89};
819	static const G4double nnFrac[10] = {0.0, 0.0, 0.2980, 0.7196, 0.9812, 0.8363, 0.5602, 0.9601, 0.9901, 1.0};
820
821	static G4CascadeInterpolator<10> interp(nnke); // Only need one!
822	pA = interp.interpolate(ekin, nnA);
823	pC = interp.interpolate(ekin, nnC);
824	pCos = interp.interpolate(ekin, nnCos);
825	pFrac = interp.interpolate(ekin, nnFrac);
826	} else if (is == 5 \|\| is == 6 \|\| is == 13 \|\| is == 34 \|\| is == 22 \|\|
827	is == 15) {
828	// pi-p, pi+n, k-p, k0bn, k+n, or k0p
829	if (verboseLevel > 3) G4cout << " meson-nucleon elastic (2)" << G4endl;
830
831	static const G4double nnke[10] = {0.0, 0.062, 0.12, 0.217, 0.533, 0.873, 1.34, 2.86, 5.86, 10.0};
832	static const G4double nnA[10] = {0.0, 27.08, 19.32, 9.92, 7.74, 9.86, 5.51, 7.25, 7.23, 7.3};
833	static const G4double nnC[10] = {0.0, 0.0, -19.49, -15.78, -9.78, -2.74, -1.16, 2.31, 2.96, 7.3};
834	static const G4double nnCos[10] = {-1.0, -1.0, -0.235, -0.259, -0.276, 0.336, 0.250, 0.732, 0.875, 0.9};
835	static const G4double nnFrac[10] = {1.0, 1.0, 0.6918, 0.6419, 0.7821, 0.6542, 0.8382, 0.9722, 0.9784, 1.0};
836
837	static G4CascadeInterpolator<10> interp(nnke); // Only need one!
838	pA = interp.interpolate(ekin, nnA);
839	pC = interp.interpolate(ekin, nnC);
840	pCos = interp.interpolate(ekin, nnCos);
841	pFrac = interp.interpolate(ekin, nnFrac);
842	} else {
843	G4cout << " G4ElementaryParticleCollider::sampleCMmomentumFor2to2: interaction not recognized "
844	<< G4endl;
845	}
846
847	// Use parameters determined above to get polar angle
848	G4double ct = sampleCMcosFor2to2(pscm, pFrac, pA, pC, pCos);
849
850	return generateWithFixedTheta(ct, pscm);
851	}
852
853
854	G4double
855	G4ElementaryParticleCollider::sampleCMcosFor2to2(G4double pscm, G4double pFrac,
856	G4double pA, G4double pC,
857	G4double pCos) const
858	{
859	G4double term1;
860	G4double term2;
861	G4double randScale;
862	G4double randVal;
863	G4double costheta = 1.0;
864
865	if (verboseLevel>3) {
866	G4cout << " sampleCMcosFor2to2: pscm " << pscm << " pFrac " << pFrac
867	<< " pA " << pA << " pC " << pC << " pCos " << pCos << G4endl;
868	}
869
870	if (G4UniformRand() < pFrac) {
871	// Sample small angles ( 0 < theta < theta0 )
872	term1 = 2.0pscmpscm*pA;
873	term2 = std::exp(-2.0*term1);
874	randScale = (std::exp(-term1*(1.0 - pCos)) - term2)/(1.0 - term2);
875	randVal = (1.0 - randScale)*G4UniformRand() + randScale;
876	} else {
877	// Sample large angles ( theta0 < theta < 180 )
878	term1 = 2.0pscmpscm*pC;
879	term2 = std::exp(-2.0*term1);
880	randScale = (std::exp(-term1*(1.0 - pCos)) - term2)/(1.0 - term2);
881	randVal = randScale*G4UniformRand();
882	}
883
884	costheta = 1.0 + std::log(randVal*(1.0 - term2) + term2)/term1;
885
886	if (verboseLevel>3) {
887	G4cout << " term1 " << term1 << " term2 " << term2 << " randVal "
888	<< randVal << " => costheta " << costheta << G4endl;
889	}
890
891	return costheta;
892	}
893
894
895	void
896	G4ElementaryParticleCollider::generateSCMpionAbsorption(G4double etot_scm,
897	G4InuclElementaryParticle* particle1,
898	G4InuclElementaryParticle* particle2) {
899	if (verboseLevel > 3) G4cout << " >>> G4ElementaryParticleCollider::generateSCMpionAbsorption"
900	<< G4endl;
901
902	// generate nucleons momenta for pion absorption
903	// the nucleon distribution assumed to be isotropic in SCM
904
905	G4InuclElementaryParticle dummy;
906
907	particles.clear(); // Initialize buffers for this event
908	particle_kinds.clear();
909
910	G4int type1 = particle1->type();
911	G4int type2 = particle2->type();
912
913	// generate kinds
914
915	if (type1 == 3) {
916
917	if (type2 == 111) { // pi+ + PP -> ?
918
919	G4cout << " pion absorption: pi+ + PP -> ? " << G4endl;
920
921	return;
922	}
923	else if(type2 == 112) { // pi+ + PN -> PP
924	particle_kinds.push_back(1);
925	particle_kinds.push_back(1);
926	}
927	else if(type2 == 122) { // pi+ + NN -> PN
928	particle_kinds.push_back(1);
929	particle_kinds.push_back(2);
930	};
931	}
932	else if(type1 == 5) {
933	if(type2 == 111) { // pi- + PP -> PN
934	particle_kinds.push_back(1);
935	particle_kinds.push_back(2);
936	}
937	else if(type2 == 112) { // pi- + PN -> NN
938	particle_kinds.push_back(2);
939	particle_kinds.push_back(2);
940	}
941	else if(type2 == 122) { // pi- + NN -> ?
942
943	G4cout << " pion absorption: pi- + NN -> ? " << G4endl;
944
945	return;
946	};
947	}
948	else if(type1 == 7) {
949	if(type2 == 111) { // pi0 + PP -> PP
950	particle_kinds.push_back(1);
951	particle_kinds.push_back(1);
952	}
953	else if(type2 == 112) { // pi0 + PN -> PN
954	particle_kinds.push_back(1);
955	particle_kinds.push_back(2);
956	}
957	else if(type2 == 122) { // pi0 + NN -> ?
958	particle_kinds.push_back(2);
959	particle_kinds.push_back(2);
960	};
961	};
962
963	G4double m1 = dummy.getParticleMass(particle_kinds[0]);
964	G4double m1sq = m1*m1;
965
966	G4double m2 = dummy.getParticleMass(particle_kinds[1]);
967	G4double m2sq = m2*m2;
968
969	G4double a = 0.5 * (etot_scm * etot_scm - m1sq - m2sq);
970
971	G4double pmod = std::sqrt((a * a - m1sq * m2sq) / (m1sq + m2sq + 2.0 * a));
972	G4LorentzVector mom1 = generateWithRandomAngles(pmod, m1);
973	G4LorentzVector mom2;
974	mom2.setVectM(-mom1.vect(), m2);
975
976	particles.push_back(G4InuclElementaryParticle(mom1, particle_kinds[0],3));
977	particles.push_back(G4InuclElementaryParticle(mom2, particle_kinds[1],3));
978
979	return;
980	}
981
982
983	// Parameter array for momentum calculation of many body final states
984	const G4double G4ElementaryParticleCollider::rmn[14][10][2] = {
985	{{0.5028, 0.6305}, {3.1442, -3.7333}, {-7.8172, 13.464}, {8.1667, -18.594},
986	{1.6208, 1.9439}, {-4.3139, -4.6268}, {12.291, 9.7879}, {-15.288, -9.6074},
987	{ 0.0, 0.0}, { 0.0, 0.0}},
988
989	{{0.9348, 2.1801}, {-10.59, 1.5163}, { 29.227, -16.38}, {-34.55, 27.944},
990	{-0.2009, -0.3464}, {1.3641, 1.1093}, {-3.403, -1.9313}, { 3.8559, 1.7064},
991	{ 0.0, 0.0}, { 0.0, 0.0}},
992
993	{{-0.0967, -1.2886}, {4.7335, -2.457}, {-14.298, 15.129}, {17.685, -23.295},
994	{ 0.0126, 0.0271}, {-0.0835, -0.1164}, { 0.186, 0.2697}, {-0.2004, -0.3185},
995	{ 0.0, 0.0}, { 0.0, 0.0}},
996
997	{{-0.025, 0.2091}, {-0.6248, 0.5228}, { 2.0282, -2.8687}, {-2.5895, 4.2688},
998	{-0.0002, -0.0007}, {0.0014, 0.0051}, {-0.0024, -0.015}, { 0.0022, 0.0196},
999	{ 0.0, 0.0}, { 0.0, 0.0}},
1000
1001	{{1.1965, 0.9336}, {-0.8289,-1.8181}, { 1.0426, 5.5157}, { -1.909,-8.5216},
1002	{ 1.2419, 1.8693}, {-4.3633, -5.5678}, { 13.743, 14.795}, {-18.592, -16.903},
1003	{ 0.0, 0.0}, { 0.0, 0.0}},
1004
1005	{{0.287, 1.7811}, {-4.9065,-8.2927}, { 16.264, 20.607}, {-19.904,-20.827},
1006	{-0.244, -0.4996}, {1.3158, 1.7874}, {-3.5691, -4.133}, { 4.3867, 3.8393},
1007	{ 0.0, 0.0}, { 0.0, 0.0}},
1008
1009	{{-0.2449, -1.5264}, { 2.9191, 6.8433}, {-9.5776, -16.067}, { 11.938, 16.845},
1010	{0.0157, 0.0462}, {-0.0826, -0.1854}, { 0.2143, 0.4531}, {-0.2585, -0.4627},
1011	{ 0.0, 0.0}, { 0.0, 0.0}},
1012
1013	{{0.0373, 0.2713}, {-0.422, -1.1944}, { 1.3883, 2.7495}, {-1.7476,-2.9045},
1014	{-0.0003, -0.0013}, {0.0014, 0.0058}, {-0.0034,-0.0146}, { 0.0039, 0.0156},
1015	{ 0.0, 0.0}, { 0.0, 0.0}},
1016
1017	{{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1018	{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1019	{ 0.1451, 0.0929},{ 0.1538, 0.1303}},
1020
1021	{{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1022	{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1023	{ 0.4652, 0.5389},{ 0.2744, 0.4071}},
1024
1025	{{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1026	{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1027	{ -0.033, -0.0545},{-0.0146, -0.0288}},
1028
1029	{{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1030	{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1031	{ 0.6296, 0.1491},{ 0.8381, 0.1802}},
1032
1033	{{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1034	{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1035	{ 0.1787, 0.385},{ 0.0086, 0.3302}},
1036
1037	{{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1038	{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1039	{-0.0026, -0.0128},{ 0.0033, -0.0094}}
1040	};
1041
1042	const G4double G4ElementaryParticleCollider::abn[4][4][4] = {
1043	{{0.0856, 0.0716, 0.1729, 0.0376}, {5.0390, 3.0960, 7.1080, 1.4331},
1044	{-13.782, -11.125, -17.961, -3.1350}, {14.661, 18.130, 16.403, 6.4864}},
1045	{{0.0543, 0.0926, -0.1450, 0.2383}, {-9.2324, -3.2186, -13.032, 1.8253},
1046	{36.397, 20.273, 41.781, 1.7648}, {-42.962, -33.245, -40.799, -16.735}},
1047	{{-0.0511, -0.0515, 0.0454, -0.1541}, {4.6003, 0.8989, 8.3515, -1.5201},
1048	{-20.534, -7.5084, -30.260, -1.5692}, {27.731, 13.188, 32.882, 17.185}},
1049	{{0.0075, 0.0058, -0.0048, 0.0250}, {-0.6253, -0.0017, -1.4095, 0.3059},
1050	{2.9159, 0.7022, 5.3505, 0.3252}, {-4.1101, -1.4854, -6.0946, -3.5277}}
1051	};

Note: See TracBrowser for help on using the repository browser.

Download in other formats: