Context Navigation

source: trunk/source/processes/hadronic/models/high_energy/src/G4HEInelastic.cc @ 962

Last change on this file since 962 was 962, checked in by garnier, 15 years ago
update processes
File size: 207.5 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	//
26	//
27	//
28
29	#include "globals.hh"
30	#include "G4ios.hh"
31
32	//
33	// G4 Process: Gheisha High Energy Collision model.
34	// This includes the high energy cascading model, the two-body-resonance model
35	// and the low energy two-body model. Not included is the low energy stuff
36	// like nuclear reactions, nuclear fission without any cascading and all
37	// processes for particles at rest.
38	//
39	// H. Fesefeldt, RWTH-Aachen, 23-October-1996
40	// Last modified: 29-July-1998
41	// HPW, fixed bug in getting pdgencoding for nuclei
42	// Hisaya, fixed HighEnergyCascading
43	// Fesefeldt, fixed bug in TuningOfHighEnergyCascading, 23 June 2000
44	// Fesefeldt, fixed next bug in TuningOfHighEnergyCascading, 14 August 2000
45	//
46	#include "G4HEInelastic.hh"
47	#include "G4HEVector.hh"
48	#include "G4ParticleDefinition.hh"
49	#include "G4DynamicParticle.hh"
50	#include "G4ParticleTable.hh"
51	#include "G4KaonZero.hh"
52	#include "G4AntiKaonZero.hh"
53	#include "G4Deuteron.hh"
54	#include "G4Triton.hh"
55	#include "G4Alpha.hh"
56
57	void G4HEInelastic::FillParticleChange(G4HEVector pv[], G4int aVecLength)
58	{
59	theParticleChange.Clear();
60	for (G4int i=0; i<aVecLength; i++)
61	{
62	G4int pdgCode = pv[i].getCode();
63	G4ParticleDefinition * aDefinition=NULL;
64	if(pdgCode == 0)
65	{
66	G4int bNumber = pv[i].getBaryonNumber();
67	if(bNumber==2) aDefinition = G4Deuteron::Deuteron();
68	if(bNumber==3) aDefinition = G4Triton::Triton();
69	if(bNumber==4) aDefinition = G4Alpha::Alpha();
70	}
71	else
72	{
73	aDefinition = G4ParticleTable::GetParticleTable()->FindParticle(pdgCode);
74	}
75	G4DynamicParticle * aParticle = new G4DynamicParticle();
76	aParticle->SetDefinition(aDefinition);
77	aParticle->SetMomentum(pv[i].getMomentum()*GeV);
78	theParticleChange.AddSecondary(aParticle);
79	G4ParticleDefinition * dummy = G4KaonZero::KaonZero();
80	dummy = G4AntiKaonZero::AntiKaonZero();
81	}
82	}
83
84	void
85	G4HEInelastic::SetParticles()
86	{
87	PionPlus.setDefinition("PionPlus");
88	PionZero.setDefinition("PionZero");
89	PionMinus.setDefinition("PionMinus");
90	KaonPlus.setDefinition("KaonPlus");
91	KaonZero.setDefinition("KaonZero");
92	AntiKaonZero.setDefinition("AntiKaonZero");
93	KaonMinus.setDefinition("KaonMinus");
94	KaonZeroShort.setDefinition("KaonZeroShort");
95	KaonZeroLong.setDefinition("KaonZeroLong");
96	Proton.setDefinition("Proton");
97	AntiProton.setDefinition("AntiProton");
98	Neutron.setDefinition("Neutron");
99	AntiNeutron.setDefinition("AntiNeutron");
100	Lambda.setDefinition("Lambda");
101	AntiLambda.setDefinition("AntiLambda");
102	SigmaPlus.setDefinition("SigmaPlus");
103	SigmaZero.setDefinition("SigmaZero");
104	SigmaMinus.setDefinition("SigmaMinus");
105	AntiSigmaPlus.setDefinition("AntiSigmaPlus");
106	AntiSigmaZero.setDefinition("AntiSigmaZero");
107	AntiSigmaMinus.setDefinition("AntiSigmaMinus");
108	XiZero.setDefinition("XiZero");
109	XiMinus.setDefinition("XiMinus");
110	AntiXiZero.setDefinition("AntiXiZero");
111	AntiXiMinus.setDefinition("AntiXiMinus");
112	OmegaMinus.setDefinition("OmegaMinus");
113	AntiOmegaMinus.setDefinition("AntiOmegaMinus");
114	Deuteron.setDefinition("Deuteron");
115	Triton.setDefinition("Triton");
116	Alpha.setDefinition("Alpha");
117	Gamma.setDefinition("Gamma");
118	return;
119	}
120
121	G4double
122	G4HEInelastic::Amin(G4double a, G4double b)
123	{
124	G4double c = a;
125	if(b < a) c = b;
126	return c;
127	}
128	G4double
129	G4HEInelastic::Amax(G4double a, G4double b)
130	{
131	G4double c = a;
132	if(b > a) c = b;
133	return c;
134	}
135	G4int
136	G4HEInelastic::Imin(G4int a, G4int b)
137	{
138	G4int c = a;
139	if(b < a) c = b;
140	return c;
141	}
142	G4int
143	G4HEInelastic::Imax(G4int a, G4int b)
144	{
145	G4int c = a;
146	if(b > a) c = b;
147	return c;
148	}
149
150
151	G4double
152	G4HEInelastic::NuclearInelasticity(G4double incidentKineticEnergy,
153	G4double atomicWeight,
154	G4double /* atomicNumber*/)
155	{
156	G4double expu = std::log(MAXFLOAT);
157	G4double expl = -expu;
158	G4double ala = std::log(atomicWeight);
159	G4double ale = std::log(incidentKineticEnergy);
160	G4double sig1 = 0.5;
161	G4double sig2 = 0.5;
162	G4double em = Amin(0.239 + 0.0408alaala, 1.);
163	G4double cinem = Amin(0.0019*std::pow(ala,3.), 0.15);
164	G4double sig = (ale > em) ? sig2 : sig1;
165	G4double corr = Amin(Amax(-std::pow(ale-em,2.)/(2.sigsig),expl), expu);
166	G4double dum1 = -(incidentKineticEnergy)*cinem;
167	G4double dum2 = std::abs(dum1);
168	G4double dum3 = std::exp(corr);
169	G4double cinema = 0.;
170	if (dum2 >= 1.) cinema = dum1*dum3;
171	else if (dum3 > 1.e-10) cinema = dum1*dum3;
172	cinema = - Amax(-incidentKineticEnergy, cinema);
173	if(verboseLevel > 1) {
174	G4cout << " NuclearInelasticity: " << ala << " " << ale << " "
175	<< em << " " << corr << " " << dum1 << " " << dum2 << " "
176	<< dum3 << " " << cinema << G4endl;
177	}
178	return cinema;
179	}
180
181	G4double
182	G4HEInelastic::NuclearExcitation(G4double incidentKineticEnergy,
183	G4double atomicWeight,
184	G4double atomicNumber,
185	G4double& excitationEnergyGPN,
186	G4double& excitationEnergyDTA)
187	{
188	G4double neutronMass = Neutron.getMass();
189	G4double electronMass = 0.000511;
190	G4double exnu = 0.;
191	excitationEnergyGPN = 0.;
192	excitationEnergyDTA = 0.;
193
194	if (atomicWeight > (neutronMass + 2.*electronMass))
195	{
196	G4int magic = ((G4int)(atomicNumber+0.1) == 82) ? 1 : 0;
197	G4double ekin = Amin(Amax(incidentKineticEnergy, 0.1), 4.);
198	G4double cfa = Amax(0.35 +((0.35 - 0.05)/2.3)*std::log(ekin), 0.15);
199	exnu = 7.716cfastd::exp(-cfa);
200	G4double atno = Amin(atomicWeight, 120.);
201	cfa = ((atno - 1.)/120.) * std::exp(-(atno-1.)/120.);
202	exnu = exnu * cfa;
203	G4double fpdiv = Amax(1.-0.25ekinekin, 0.5);
204	G4double gfa = 2.*((atomicWeight-1.)/70.)
205	* std::exp(-(atomicWeight-1.)/70.);
206
207	excitationEnergyGPN = exnu * fpdiv;
208	excitationEnergyDTA = exnu - excitationEnergyGPN;
209
210	G4double ran1 = 0., ran2 = 0.;
211	if (!magic)
212	{ ran1 = normal();
213	ran2 = normal();
214	}
215	excitationEnergyGPN = Amax(excitationEnergyGPN(1.+ran1gfa),0.);
216	excitationEnergyDTA = Amax(excitationEnergyDTA(1.+ran2gfa),0.);
217	exnu = excitationEnergyGPN + excitationEnergyDTA;
218	if(verboseLevel > 1) {
219	G4cout << " NuclearExcitation: " << magic << " " << ekin
220	<< " " << cfa << " " << atno << " " << fpdiv << " "
221	<< gfa << " " << excitationEnergyGPN
222	<< " " << excitationEnergyDTA << G4endl;
223	}
224
225	while (exnu >= incidentKineticEnergy)
226	{
227	excitationEnergyGPN = (1. - 0.5normal());
228	excitationEnergyDTA = (1. - 0.5normal());
229	exnu = excitationEnergyGPN + excitationEnergyDTA;
230	}
231	}
232	return exnu;
233	}
234
235	G4double
236	G4HEInelastic::pmltpc(G4int np, G4int nm, G4int nz, G4int n,
237	G4double b, G4double c)
238	{
239	G4double expxu = std::log(MAXFLOAT);
240	G4double expxl = -expxu;
241	G4int i;
242	G4double npf = 0.0, nmf = 0.0, nzf = 0.0;
243	for(i=2;i<=np;i++) npf += std::log((G4double)i);
244	for(i=2;i<=nm;i++) nmf += std::log((G4double)i);
245	for(i=2;i<=nz;i++) nzf += std::log((G4double)i);
246	G4double r = Amin(expxu,Amax(expxl,-(np-nm+nz+b)(np-nm+nz+b)/(2ccn*n)-npf-nmf-nzf));
247	return std::exp(r);
248	}
249
250
251	G4int G4HEInelastic::Factorial(G4int n)
252	{
253	G4int result = 1;
254	if (n < 0) G4Exception("G4HEInelastic::Factorial()", "601",
255	FatalException, "Negative factorial argument");
256	while (n > 1) result *= n--;
257	return result;
258	}
259
260
261	G4double G4HEInelastic::normal()
262	{
263	G4double ran = -6.0;
264	for(G4int i=0; i<12; i++) ran += G4UniformRand();
265	return ran;
266	}
267
268	G4int G4HEInelastic::Poisson( G4double x )
269	{
270	G4int i, iran = 0;
271	G4double ran;
272	if ( x > 9.9 )
273	{
274	iran = G4int( Amax( 0.0, x + normal() * std::sqrt( x ) ) );
275	}
276	else
277	{
278	G4int mm = G4int( 5.0 * x );
279	if ( mm <= 0 )
280	{
281	G4double p1 = x * std::exp( -x );
282	G4double p2 = x * p1/2.;
283	G4double p3 = x * p2/3.;
284	ran = G4UniformRand();
285	if ( ran < p3 ) iran = 3;
286	else if ( ran < p2 ) iran = 2;
287	else if ( ran < p1 ) iran = 1;
288	}
289	else
290	{ G4double r = std::exp( -x );
291	ran = G4UniformRand();
292	if (ran > r)
293	{
294	G4double rrr;
295	G4double rr = r;
296	for (i=1; i <= mm; i++)
297	{
298	iran++;
299	if ( i > 5 ) rrr = std::exp(istd::log(x)-(i+0.5)std::log((G4double)i)+i-0.9189385);
300	else rrr = std::pow(x,i)*Factorial(i);
301	rr += r * rrr;
302	if (ran <= rr) break;
303	}
304	}
305	}
306	}
307	return iran;
308	}
309	G4double
310	G4HEInelastic::GammaRand( G4double avalue )
311	{
312	G4double ga = avalue -1.;
313	G4double la = std::sqrt(2.*avalue - 1.);
314	G4double ep = 1.570796327 + std::atan(ga/la);
315	G4double ro = 1.570796327 - ep;
316	G4double y = 1.;
317	G4double xtrial;
318	repeat:
319	xtrial = ga + la * std::tan(ep*G4UniformRand() + ro);
320	if(xtrial == 0.) goto repeat;
321	y = std::log(1.+sqr((xtrial-ga)/la))+ga*std::log(xtrial/ga)-xtrial+ga;
322	if(std::log(G4UniformRand()) > y) goto repeat;
323	return xtrial;
324	}
325	G4double
326	G4HEInelastic::Erlang( G4int mvalue )
327	{
328	G4double ran = G4UniformRand();
329	G4double xtrial = 0.62666*std::log((1.+ran)/(1.-ran));
330	if(G4UniformRand()<0.5) xtrial = -xtrial;
331	return mvalue+xtrial*std::sqrt(G4double(mvalue));
332	}
333
334	void
335	G4HEInelastic::StrangeParticlePairProduction(
336	const G4double availableEnergy,
337	const G4double centerOfMassEnergy,
338	G4HEVector pv[],
339	G4int &vecLen,
340	G4HEVector incidentParticle,
341	G4HEVector targetParticle )
342
343	// Choose charge combinations K+ K-, K+ K0, K0 K0, K0 K-,
344	// K+ Y0, K0 Y+, K0 Y-
345	// For antibaryon induced reactions half of the cross sections KB YB
346	// pairs are produced. Charge is not conserved, no experimental data
347	// available for exclusive reactions, therefore some average behavior
348	// assumed. The ratio L/SIGMA is taken as 3:1 (from experimental low
349	// energy data)
350
351	{
352	static G4double avrs[] = {3.,4.,5.,6.,7.,8.,9.,10.,20.,30.,40.,50.};
353	static G4double avkkb[] = {0.0015,0.0050,0.0120,0.0285,0.0525,0.0750,0.0975,
354	0.1230,0.2800,0.3980,0.4950,0.5730};
355	static G4double kkb[] = {0.2500,0.3750,0.5000,0.5625,0.6250,0.6875,0.7500,
356	0.8750,1.0000};
357	static G4double ky[] = {0.2000,0.3000,0.4000,0.5500,0.6250,0.7000,0.8000,
358	0.8500,0.9000,0.9500,0.9750,1.0000};
359	static G4int ipakkb[] = {10,13,10,11,10,12,11,11,11,12,12,11,12,12,
360	11,13,12,13};
361	static G4int ipaky[] = {18,10,18,11,18,12,20,10,20,11,20,12,21,10,
362	21,11,21,12,22,10,22,11,22,12};
363	static G4int ipakyb[] = {19,13,19,12,19,11,23,13,23,12,23,11,24,13,
364	24,12,24,11,25,13,25,12,25,11};
365	static G4double avky[] = {0.0050,0.0300,0.0640,0.0950,0.1150,0.1300,0.1450,
366	0.1550,0.2000,0.2050,0.2100,0.2120};
367	static G4double avnnb[] ={0.00001,0.0001,0.0006,0.0025,0.0100,0.0200,0.0400,
368	0.0500,0.1200,0.1500,0.1800,0.2000};
369
370	G4int i, ibin, i3, i4; // misc. local variables
371	G4double avk, avy, avn, ran;
372
373	G4double protonMass = Proton.getMass();
374	G4double sigmaMinusMass = SigmaMinus.getMass();
375	G4int antiprotonCode = AntiProton.getCode();
376	G4int antineutronCode = AntiNeutron.getCode();
377	G4int antilambdaCode = AntiLambda.getCode();
378
379	G4double incidentMass = incidentParticle.getMass();
380	G4int incidentCode = incidentParticle.getCode();
381
382	G4double targetMass = targetParticle.getMass();
383
384	// protection against annihilation processes like pbar p -> pi pi.
385
386	if (vecLen <= 2) return;
387
388	// determine the center of mass energy bin
389
390	i = 1;
391	while ( (i<12) && (centerOfMassEnergy > avrs[i]) )i++;
392	if ( i == 12 ) ibin = 11;
393	else ibin = i;
394
395	// the fortran code chooses a random replacement of produced kaons
396	// but does not take into account charge conservation
397
398	if( vecLen == 3 ) { // we know that vecLen > 2
399	i3 = 2;
400	i4 = 3; // note that we will be adding a new
401	} // secondary particle in this case only
402	else
403	{ // otherwise 2 <= i3,i4 <= vecLen
404	i4 = i3 = 2 + G4int( (vecLen-2)*G4UniformRand() );
405	while ( i3 == i4 ) i4 = 2 + G4int( (vecLen-2)*G4UniformRand() );
406	}
407
408	// use linear interpolation or extrapolation by y=centerofmassEnergy*x+b
409
410	avk = (std::log(avkkb[ibin])-std::log(avkkb[ibin-1]))*(centerOfMassEnergy-avrs[ibin-1])
411	/(avrs[ibin]-avrs[ibin-1]) + std::log(avkkb[ibin-1]);
412	avk = std::exp(avk);
413
414	avy = (std::log(avky[ibin])-std::log(avky[ibin-1]))*(centerOfMassEnergy-avrs[ibin-1])
415	/(avrs[ibin]-avrs[ibin-1]) + std::log(avky[ibin-1]);
416	avy = std::exp(avy);
417
418	avn = (std::log(avnnb[ibin])-std::log(avnnb[ibin-1]))*(centerOfMassEnergy-avrs[ibin-1])
419	/(avrs[ibin]-avrs[ibin-1]) + std::log(avnnb[ibin-1]);
420	avn = std::exp(avn);
421
422	if ( avk+avy+avn <= 0.0 ) return;
423
424	if ( incidentMass < protonMass ) avy /= 2.0;
425	avy += avk+avn;
426	avk += avn;
427
428	ran = G4UniformRand();
429	if ( ran < avn )
430	{ // p pbar && n nbar production
431	if ( availableEnergy < 2.0) return;
432	if ( vecLen == 3 )
433	{ // add a new secondary
434	if ( G4UniformRand() < 0.5 )
435	{
436	pv[i3] = Neutron;;
437	pv[vecLen++] = AntiNeutron;
438	}
439	else
440	{
441	pv[i3] = Proton;
442	pv[vecLen++] = AntiProton;
443	}
444	}
445	else
446	{ // replace two secondaries
447	if ( G4UniformRand() < 0.5 )
448	{
449	pv[i3] = Neutron;
450	pv[i4] = AntiNeutron;
451	}
452	else
453	{
454	pv[i3] = Proton;
455	pv[i4] = AntiProton;
456	}
457	}
458	}
459	else if ( ran < avk )
460	{ // K Kbar production
461	if ( availableEnergy < 1.0) return;
462	G4double ran1 = G4UniformRand();
463	i = 0;
464	while( (i<9) && (ran1>kkb[i]) )i++;
465	if ( i == 9 ) return;
466
467	// ipakkb[] = { 10,13, 10,11, 10,12, 11, 11, 11,12, 12,11, 12,12, 11,13, 12,13 };
468	// charge K+ K- K+ K0S K+ K0L K0S K0S K0S K0L K0LK0S K0LK0L K0S K- K0LK-
469
470	switch( ipakkb[i*2] )
471	{
472	case 10: pv[i3] = KaonPlus; break;
473	case 11: pv[i3] = KaonZeroShort;break;
474	case 12: pv[i3] = KaonZeroLong; break;
475	case 13: pv[i3] = KaonMinus; break;
476	}
477
478	if( vecLen == 2 )
479	{ // add a secondary
480	switch( ipakkb[i*2+1] )
481	{
482	case 10: pv[vecLen++] = KaonPlus; break;
483	case 11: pv[vecLen++] = KaonZeroShort;break;
484	case 12: pv[vecLen++] = KaonZeroLong; break;
485	case 13: pv[vecLen++] = KaonMinus; break;
486	}
487	}
488	else
489	{ // replace
490	switch( ipakkb[i*2+1] )
491	{
492	case 10: pv[i4] = KaonPlus; break;
493	case 11: pv[i4] = KaonZeroShort;break;
494	case 12: pv[i4] = KaonZeroLong; break;
495	case 13: pv[i4] = KaonMinus; break;
496	}
497	}
498	}
499	else if ( ran < avy )
500	{ // Lambda K && Sigma K
501	if( availableEnergy < 1.6) return;
502	G4double ran1 = G4UniformRand();
503	i = 0;
504	while( (i<12) && (ran1>ky[i]) )i++;
505	if ( i == 12 ) return;
506	if ( (incidentMass<protonMass) \|\| (G4UniformRand()<0.5) )
507	{
508
509	// ipaky[] = { 18,10, 18,11, 18,12, 20,10, 20,11, 20,12,
510	// L0 K+ L0 K0S L0 K0L S+ K+ S+ K0S S+ K0L
511	//
512	// 21,10, 21,11, 21,12, 22,10, 22,11, 22,12 }
513	// S0 K+ S0 K0S S0 K0L S- K+ S- K0S S- K0L
514
515	switch( ipaky[i*2] )
516	{
517	case 18: pv[1] = Lambda; break;
518	case 20: pv[1] = SigmaPlus; break;
519	case 21: pv[1] = SigmaZero; break;
520	case 22: pv[1] = SigmaMinus;break;
521	}
522	switch( ipaky[i*2+1] )
523	{
524	case 10: pv[i3] = KaonPlus; break;
525	case 11: pv[i3] = KaonZeroShort;break;
526	case 12: pv[i3] = KaonZeroLong; break;
527	}
528	}
529	else
530	{ // Lbar K && Sigmabar K production
531
532	// ipakyb[] = { 19,13, 19,12, 19,11, 23,13, 23,12, 23,11,
533	// Lb K- Lb K0L Lb K0S S+b K- S+b K0L S+b K0S
534	// 24,13, 24,12, 24,11, 25,13, 25,12, 25,11 };
535	// S0b K- S0BK0L S0BK0S S-BK- S-B K0L S-BK0S
536
537	if( (incidentCode==antiprotonCode) \|\| (incidentCode==antineutronCode) \|\|
538	(incidentCode==antilambdaCode) \|\| (incidentMass>sigmaMinusMass) )
539	{
540	switch( ipakyb[i*2] )
541	{
542	case 19:pv[0] = AntiLambda; break;
543	case 23:pv[0] = AntiSigmaPlus; break;
544	case 24:pv[0] = AntiSigmaZero; break;
545	case 25:pv[0] = AntiSigmaMinus;break;
546	}
547	switch( ipakyb[i*2+1] )
548	{
549	case 11:pv[i3] = KaonZeroShort;break;
550	case 12:pv[i3] = KaonZeroLong; break;
551	case 13:pv[i3] = KaonMinus; break;
552	}
553	}
554	else
555	{
556	switch( ipaky[i*2] )
557	{
558	case 18:pv[0] = Lambda; break;
559	case 20:pv[0] = SigmaPlus; break;
560	case 21:pv[0] = SigmaZero; break;
561	case 22:pv[0] = SigmaMinus;break;
562	}
563	switch( ipaky[i*2+1] )
564	{
565	case 10: pv[i3] = KaonPlus; break;
566	case 11: pv[i3] = KaonZeroShort;break;
567	case 12: pv[i3] = KaonZeroLong; break;
568	}
569	}
570	}
571	}
572	else
573	return;
574
575	// check the available energy
576	// if there is not enough energy for kkb/ky pair production
577	// then reduce the number of secondary particles
578	// NOTE:
579	// the number of secondaries may have been changed
580	// the incident and/or target particles may have changed
581	// charge conservation is ignored (as well as strangness conservation)
582
583	incidentMass = incidentParticle.getMass();
584	targetMass = targetParticle.getMass();
585
586	G4double energyCheck = centerOfMassEnergy-(incidentMass+targetMass);
587	if (verboseLevel > 1) G4cout << "Particles produced: " ;
588
589	for ( i=0; i < vecLen; i++ )
590	{
591	energyCheck -= pv[i].getMass();
592	if (verboseLevel > 1) G4cout << pv[i].getCode() << " " ;
593	if( energyCheck < 0.0 )
594	{
595	if( i > 0 ) vecLen = --i; // chop off the secondary list
596	return;
597	}
598	}
599	if (verboseLevel > 1) G4cout << G4endl;
600	return;
601	}
602
603	void
604	G4HEInelastic::HighEnergyCascading(G4bool &successful,
605	G4HEVector pv[],
606	G4int &vecLen,
607	G4double &excitationEnergyGNP,
608	G4double &excitationEnergyDTA,
609	G4HEVector incidentParticle,
610	G4HEVector targetParticle,
611	G4double atomicWeight,
612	G4double atomicNumber)
613	{
614	//
615	// The multiplicity of particles produced in the first interaction has been
616	// calculated in one of the FirstIntInNuc.... routines. The nuclear
617	// cascading particles are parameterized from experimental data.
618	// A simple single variable description E D3S/DP3= F(Q) with
619	// Q^2 = (M*X)^2 + PT^2 is used. Final state kinematics are produced
620	// by an FF-type iterative cascade method.
621	// Nuclear evaporation particles are added at the end of the routine.
622
623	// All quantities in the G4HEVector Array pv are in GeV- units.
624	// The method is a copy of MediumEnergyCascading with some special tuning
625	// for high energy interactions.
626
627
628	G4int protonCode = Proton.getCode();
629	G4double protonMass = Proton.getMass();
630	G4int neutronCode = Neutron.getCode();
631	G4double neutronMass = Neutron.getMass();
632	G4double kaonPlusMass = KaonPlus.getMass();
633	G4int kaonPlusCode = KaonPlus.getCode();
634	G4int kaonMinusCode = KaonMinus.getCode();
635	G4int kaonZeroSCode = KaonZeroShort.getCode();
636	G4int kaonZeroLCode = KaonZeroLong.getCode();
637	G4int kaonZeroCode = KaonZero.getCode();
638	G4int antiKaonZeroCode = AntiKaonZero.getCode();
639	G4int pionPlusCode = PionPlus.getCode();
640	G4int pionZeroCode = PionZero.getCode();
641	G4int pionMinusCode = PionMinus.getCode();
642	G4String mesonType = PionPlus.getType();
643	G4String baryonType = Proton.getType();
644	G4String antiBaryonType= AntiProton.getType();
645
646	G4double targetMass = targetParticle.getMass();
647
648	G4int incidentCode = incidentParticle.getCode();
649	G4double incidentMass = incidentParticle.getMass();
650	G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
651	G4double incidentEnergy = incidentParticle.getEnergy();
652	G4double incidentKineticEnergy = incidentParticle.getKineticEnergy();
653	G4String incidentType = incidentParticle.getType();
654	// G4double incidentTOF = incidentParticle.getTOF();
655	G4double incidentTOF = 0.;
656
657	// some local variables
658
659	G4int i, j, l;
660
661	if (verboseLevel > 1)
662	G4cout << " G4HEInelastic::HighEnergyCascading " << G4endl;
663	successful = false;
664	if(incidentTotalMomentum < 25. + G4UniformRand()*25.) return;
665
666	// define annihilation channels.
667
668	G4bool annihilation = false;
669	if (incidentCode < 0 && incidentType == antiBaryonType &&
670	pv[0].getType() != antiBaryonType &&
671	pv[1].getType() != antiBaryonType )
672	{
673	annihilation = true;
674	}
675
676
677
678	G4double twsup[] = { 1., 1., 0.7, 0.5, 0.3, 0.2, 0.1, 0.0 };
679
680	if( annihilation ) goto start;
681	if( vecLen >= 8) goto start;
682	if( incidentKineticEnergy < 1.) return;
683	if( ( incidentCode == kaonPlusCode \|\| incidentCode == kaonMinusCode
684	\|\| incidentCode == kaonZeroCode \|\| incidentCode == antiKaonZeroCode
685	\|\| incidentCode == kaonZeroSCode \|\| incidentCode == kaonZeroLCode )
686	&& ( G4UniformRand() < 0.5) ) goto start;
687	if( G4UniformRand() > twsup[vecLen-1]) goto start;
688	if( incidentKineticEnergy > (G4UniformRand()*200 + 50.) ) goto start;
689	return;
690
691	start:
692
693	if (annihilation)
694	{ // do some corrections of incident particle kinematic
695	G4double ekcor = Amax( 1., 1./incidentKineticEnergy);
696	incidentKineticEnergy = 2targetMass + incidentKineticEnergy(1.+ekcor/atomicWeight);
697	G4double excitation = NuclearExcitation(incidentKineticEnergy,
698	atomicWeight,
699	atomicNumber,
700	excitationEnergyGNP,
701	excitationEnergyDTA);
702	incidentKineticEnergy -= excitation;
703	if (incidentKineticEnergy < excitationEnergyDTA) incidentKineticEnergy = 0.;
704	incidentEnergy = incidentKineticEnergy + incidentMass;
705	incidentTotalMomentum =
706	std::sqrt( Amax(0., incidentEnergyincidentEnergy - incidentMassincidentMass));
707	}
708
709	G4HEVector pTemp;
710	for (i = 2; i<vecLen; i++)
711	{
712	j = Imin( vecLen-1, (G4int)(2. + G4UniformRand()*(vecLen - 2)));
713	pTemp = pv[j];
714	pv[j] = pv[i];
715	pv[i] = pTemp;
716	}
717	// randomize the first two leading particles
718	// for kaon induced reactions only
719	// (need from experimental data)
720
721	if( (incidentCode==kaonPlusCode \|\| incidentCode==kaonMinusCode \|\|
722	incidentCode==kaonZeroCode \|\| incidentCode==antiKaonZeroCode \|\|
723	incidentCode==kaonZeroSCode \|\| incidentCode==kaonZeroLCode)
724	&& (G4UniformRand() > 0.9) )
725	{
726	pTemp = pv[1];
727	pv[1] = pv[0];
728	pv[0] = pTemp;
729	}
730	// mark leading particles for incident strange particles
731	// and antibaryons, for all other we assume that the first
732	// and second particle are the leading particles.
733	// We need this later for kinematic aspects of strangeness
734	// conservation.
735
736	G4int lead = 0;
737	G4HEVector leadParticle;
738	if( (incidentMass >= kaonPlusMass-0.05) && (incidentCode != protonCode)
739	&& (incidentCode != neutronCode) )
740	{
741	G4double pMass = pv[0].getMass();
742	G4int pCode = pv[0].getCode();
743	if( (pMass >= kaonPlusMass-0.05) && (pCode != protonCode)
744	&& (pCode != neutronCode) )
745	{
746	lead = pCode;
747	leadParticle = pv[0];
748	}
749	else
750	{
751	pMass = pv[1].getMass();
752	pCode = pv[1].getCode();
753	if( (pMass >= kaonPlusMass-0.05) && (pCode != protonCode)
754	&& (pCode != neutronCode) )
755	{
756	lead = pCode;
757	leadParticle = pv[1];
758	}
759	}
760	}
761
762	// Distribute particles in forward and backward hemispheres in center
763	// of mass system. Incident particle goes in forward hemisphere.
764
765	G4HEVector pvI = incidentParticle; // for the incident particle
766	pvI.setSide( 1 );
767
768	G4HEVector pvT = targetParticle; // for the target particle
769	pvT.setMomentumAndUpdate( 0.0, 0.0, 0.0 );
770	pvT.setSide( -1 );
771	pvT.setTOF( -1.);
772
773
774	G4double centerOfMassEnergy = std::sqrt( sqr(incidentMass)+sqr(targetMass)
775	+2.0targetMassincidentEnergy );
776	// G4double availableEnergy = centerOfMassEnergy - ( targetMass + incidentMass );
777
778	G4double tavai1 = centerOfMassEnergy/2.0 - incidentMass;
779	G4double tavai2 = centerOfMassEnergy/2.0 - targetMass;
780
781	// define G4HEVector- array for kinematic manipulations,
782	// with a one by one correspondence to the pv-Array.
783
784	G4int ntb = 1;
785	for( i=0; i < vecLen; i++ )
786	{
787	if (i == 0) pv[i].setSide( 1 );
788	else if (i == 1) pv[i].setSide( -1 );
789	else
790	{ if( G4UniformRand() < 0.5 )
791	{
792	pv[i].setSide( -1 );
793	ntb++;
794	}
795	else
796	pv[i].setSide( 1 );
797	}
798	pv[i].setTOF( incidentTOF);
799	}
800	G4double tb = 2. * ntb;
801	if (centerOfMassEnergy < (2. + G4UniformRand()))
802	tb = (2. * ntb + vecLen)/2.;
803
804	if (verboseLevel > 1)
805	{ G4cout << " pv Vector after Randomization " << vecLen << G4endl;
806	pvI.Print(-1);
807	pvT.Print(-1);
808	for (i=0; i < vecLen ; i++) pv[i].Print(i);
809	}
810
811	// Add particles from intranuclear cascade
812	// nuclearCascadeCount = number of new secondaries produced by nuclear
813	// cascading.
814	// extraCount = number of nucleons within these new secondaries
815
816	G4double s, xtarg, ran;
817	s = centerOfMassEnergy*centerOfMassEnergy;
818	G4double afc;
819	afc = Amin(0.5, 0.312 + 0.200 * std::log(std::log(s))+ std::pow(s,1.5)/6000.0);
820	xtarg = Amax(0.01, afc * (std::pow(atomicWeight, 0.33) - 1.0) * tb);
821	G4int nstran = Poisson( 0.03*xtarg);
822	G4int momentumBin = 0;
823	G4double nucsup[] = { 1.00, 0.7, 0.5, 0.4, 0.5, 0.5 };
824	G4double psup[] = { 3., 6., 20., 50., 100., 1000. };
825	while( (momentumBin < 6) && (incidentTotalMomentum > psup[momentumBin])) momentumBin++;
826	momentumBin = Imin(5, momentumBin);
827	G4double xpnhmf = Amax(0.01,xtarg*nucsup[momentumBin]);
828	G4double xshhmf = Amax(0.01,xtarg - xpnhmf);
829	G4double rshhmf = 0.25*xshhmf;
830	G4double rpnhmf = 0.81*xpnhmf;
831	G4double xhmf=0;
832	if(verboseLevel > 1)
833	G4cout << "xtarg= " << xtarg << " xpnhmf = " << xpnhmf << G4endl;
834
835	G4int nshhmf, npnhmf;
836	if (rshhmf > 1.1)
837	{
838	rshhmf = xshhmf/(rshhmf-1.);
839	if (rshhmf <= 20.)
840	xhmf = GammaRand( rshhmf );
841	else
842	xhmf = Erlang( G4int(rshhmf+0.5) );
843	xshhmf *= xhmf/rshhmf;
844	}
845	nshhmf = Poisson( xshhmf );
846	if(verboseLevel > 1)
847	G4cout << "xshhmf = " << xshhmf << " xhmf = " << xhmf
848	<< " rshhmf = " << rshhmf << G4endl;
849
850	if (rpnhmf > 1.1)
851	{
852	rpnhmf = xpnhmf/(rpnhmf -1.);
853	if (rpnhmf <= 20.)
854	xhmf = GammaRand( rpnhmf );
855	else
856	xhmf = Erlang( G4int(rpnhmf+0.5) );
857	xpnhmf *= xhmf/rpnhmf;
858	}
859	npnhmf = Poisson( xpnhmf );
860	if(verboseLevel > 1)
861	G4cout << "nshhmf = " << nshhmf << " npnhmf = " << npnhmf
862	<< " nstran = " << nstran << G4endl;
863
864	G4int ntarg = nshhmf + npnhmf + nstran;
865
866	G4int targ = 0;
867
868	while (npnhmf > 0)
869	{
870	if ( G4UniformRand() > (1. - atomicNumber/atomicWeight))
871	pv[vecLen] = Proton;
872	else
873	pv[vecLen] = Neutron;
874	targ++;
875	pv[vecLen].setSide( -2 );
876	pv[vecLen].setFlag( true );
877	pv[vecLen].setTOF( incidentTOF );
878	vecLen++;
879	npnhmf--;
880	}
881	while (nstran > 0)
882	{
883	ran = G4UniformRand();
884	if (ran < 0.14) pv[vecLen] = Lambda;
885	else if (ran < 0.20) pv[vecLen] = SigmaZero;
886	else if (ran < 0.43) pv[vecLen] = KaonPlus;
887	else if (ran < 0.66) pv[vecLen] = KaonZero;
888	else if (ran < 0.89) pv[vecLen] = AntiKaonZero;
889	else pv[vecLen] = KaonMinus;
890	if (G4UniformRand() > 0.2)
891	{
892	pv[vecLen].setSide( -2 );
893	pv[vecLen].setFlag( true );
894	}
895	else
896	{
897	pv[vecLen].setSide( 1 );
898	pv[vecLen].setFlag( false );
899	ntarg--;
900	}
901	pv[vecLen].setTOF( incidentTOF );
902	vecLen++;
903	nstran--;
904	}
905	while (nshhmf > 0)
906	{
907	ran = G4UniformRand();
908	if( ran < 0.33333 )
909	pv[vecLen] = PionPlus;
910	else if( ran < 0.66667 )
911	pv[vecLen] = PionZero;
912	else
913	pv[vecLen] = PionMinus;
914	if (G4UniformRand() > 0.2)
915	{
916	pv[vecLen].setSide( -2 ); // backward cascade particles
917	pv[vecLen].setFlag( true ); // true is the same as IPA(i)<0
918	}
919	else
920	{
921	pv[vecLen].setSide( 1 );
922	pv[vecLen].setFlag( false );
923	ntarg--;
924	}
925	pv[vecLen].setTOF( incidentTOF );
926	vecLen++;
927	nshhmf--;
928	}
929
930	// assume conservation of kinetic energy
931	// in forward & backward hemispheres
932
933	G4int is, iskip, iavai1;
934	if(vecLen <= 1) return;
935
936	tavai1 = centerOfMassEnergy/2.;
937	iavai1 = 0;
938
939	for (i = 0; i < vecLen; i++)
940	{
941	if (pv[i].getSide() > 0)
942	{
943	tavai1 -= pv[i].getMass();
944	iavai1++;
945	}
946	}
947	if ( iavai1 == 0) return;
948
949	while( tavai1 <= 0.0 )
950	{ // must eliminate a particle from the forward side
951	iskip = G4int(G4UniformRand()*iavai1) + 1;
952	is = 0;
953	for( i=vecLen-1; i>=0; i-- )
954	{
955	if( pv[i].getSide() > 0 )
956	{
957	if (++is == iskip)
958	{
959	tavai1 += pv[i].getMass();
960	iavai1--;
961	if ( i != vecLen-1)
962	{
963	for( j=i; j<vecLen; j++ )
964	{
965	pv[j] = pv[j+1];
966	}
967	}
968	if( --vecLen == 0 ) return; // all the secondaries except of the
969	break; // --+
970	} // \|
971	} // v
972	} // break goes down to here
973	} // to the end of the for- loop.
974
975
976	tavai2 = (targ+1)*centerOfMassEnergy/2.;
977	G4int iavai2 = 0;
978
979	for (i = 0; i < vecLen; i++)
980	{
981	if (pv[i].getSide() < 0)
982	{
983	tavai2 -= pv[i].getMass();
984	iavai2++;
985	}
986	}
987	if (iavai2 == 0) return;
988
989	while( tavai2 <= 0.0 )
990	{ // must eliminate a particle from the backward side
991	iskip = G4int(G4UniformRand()*iavai2) + 1;
992	is = 0;
993	for( i = vecLen-1; i >= 0; i-- )
994	{
995	if( pv[i].getSide() < 0 )
996	{
997	if( ++is == iskip )
998	{
999	tavai2 += pv[i].getMass();
1000	iavai2--;
1001	if (pv[i].getSide() == -2) ntarg--;
1002	if (i != vecLen-1)
1003	{
1004	for( j=i; j<vecLen; j++)
1005	{
1006	pv[j] = pv[j+1];
1007	}
1008	}
1009	if (--vecLen == 0) return;
1010	break;
1011	}
1012	}
1013	}
1014	}
1015
1016	if (verboseLevel > 1)
1017	{ G4cout << " pv Vector after Energy checks "
1018	<< vecLen << " " << tavai1 << " " << iavai1 << " " << tavai2
1019	<< " " << iavai2 << " " << ntarg << G4endl;
1020	pvI.Print(-1);
1021	pvT.Print(-1);
1022	for (i=0; i < vecLen ; i++) pv[i].Print(i);
1023	}
1024
1025	// define some vectors for Lorentz transformations
1026
1027	G4HEVector* pvmx = new G4HEVector [10];
1028
1029	pvmx[0].setMass( incidentMass );
1030	pvmx[0].setMomentumAndUpdate( 0.0, 0.0, incidentTotalMomentum );
1031	pvmx[1].setMass( protonMass);
1032	pvmx[1].setMomentumAndUpdate( 0.0, 0.0, 0.0 );
1033	pvmx[3].setMass( protonMass*(1+targ));
1034	pvmx[3].setMomentumAndUpdate( 0.0, 0.0, 0.0 );
1035	pvmx[4].setZero();
1036	pvmx[5].setZero();
1037	pvmx[7].setZero();
1038	pvmx[8].setZero();
1039	pvmx[8].setMomentum( 1.0, 0.0 );
1040	pvmx[2].Add( pvmx[0], pvmx[1] );
1041	pvmx[3].Add( pvmx[3], pvmx[0] );
1042	pvmx[0].Lor( pvmx[0], pvmx[2] );
1043	pvmx[1].Lor( pvmx[1], pvmx[2] );
1044
1045	if (verboseLevel > 1)
1046	{ G4cout << " General Vectors after Definition " << G4endl;
1047	for (i=0; i<10; i++) pvmx[i].Print(i);
1048	}
1049
1050	// Main loop for 4-momentum generation - see Pitha-report (Aachen)
1051	// for a detailed description of the method.
1052	// Process the secondary particles in reverse order.
1053
1054	G4double dndl[20];
1055	G4double binl[20];
1056	G4double pvMass(0), pvEnergy(0);
1057	G4int pvCode;
1058	G4double aspar, pt, phi, et, xval;
1059	G4double ekin = 0.;
1060	G4double ekin1 = 0.;
1061	G4double ekin2 = 0.;
1062	G4int npg = 0;
1063	G4double rmg0 = 0.;
1064	G4int targ1 = 0; // No fragmentation model for nucleons from
1065	phi = G4UniformRand()*twopi;
1066	for( i=vecLen-1; i>=0; i-- ) // the intranuclear cascade. Mark them with
1067	{ // -3 and leave the loop
1068	if( i == 1)
1069	{
1070	if ( (pv[i].getMass() > neutronMass + 0.05) && (G4UniformRand() < 0.2))
1071	{
1072	if(++npg < 19)
1073	{
1074	pv[i].setSide(-3);
1075	rmg0 += pv[i].getMass();
1076	targ++;
1077	continue;
1078	}
1079	}
1080	else if ( pv[i].getMass() > protonMass - 0.05)
1081	{
1082	if(++npg < 19)
1083	{
1084	pv[i].setSide(-3);
1085	rmg0 += pv[i].getMass();
1086	targ++;
1087	continue;
1088	}
1089	}
1090	}
1091	if( pv[i].getSide() == -2)
1092	{
1093	if ( pv[i].getName() == "Proton" \|\| pv[i].getName() == "Neutron")
1094	{
1095	if( ++npg < 19 )
1096	{
1097	pv[i].setSide( -3 );
1098	rmg0 += pv[i].getMass();
1099	targ1++;
1100	continue; // leave the for loop !!
1101	}
1102	}
1103	}
1104	// Set pt and phi values - they are changed somewhat in the
1105	// iteration loop.
1106	// Set mass parameter for lambda fragmentation model
1107
1108	G4double maspar[] = { 0.75, 0.70, 0.65, 0.60, 0.50, 0.40, 0.20, 0.10};
1109	G4double bp[] = { 4.00, 2.50, 2.20, 3.00, 3.00, 1.70, 3.50, 3.50};
1110	G4double ptex[] = { 1.70, 1.70, 1.50, 1.70, 1.40, 1.20, 1.70, 1.20};
1111
1112	// Set parameters for lambda simulation
1113	// pt is the average transverse momentum
1114	// aspar is average transverse mass
1115
1116	pvMass = pv[i].getMass();
1117	j = 2;
1118	if (pv[i].getType() == mesonType ) j = 1;
1119	if ( pv[i].getMass() < 0.4 ) j = 0;
1120	if ( i <= 1 ) j += 3;
1121	if (pv[i].getSide() <= -2) j = 6;
1122	if (j == 6 && (pv[i].getType() == baryonType \|\| pv[i].getType() == antiBaryonType)) j = 7;
1123	pt = std::sqrt(std::pow(-std::log(1.-G4UniformRand())/bp[j],ptex[j]));
1124	if(pt<0.05) pt = Amax(0.001, 0.3*G4UniformRand());
1125	aspar = maspar[j];
1126	phi = G4UniformRand()*twopi;
1127	pv[i].setMomentum( ptstd::cos(phi), ptstd::sin(phi) ); // set x- and y-momentum
1128
1129	for( j=0; j<20; j++ ) binl[j] = j/(19.*pt); // set the lambda - bins.
1130
1131	if( pv[i].getSide() > 0 )
1132	et = pvmx[0].getEnergy();
1133	else
1134	et = pvmx[1].getEnergy();
1135
1136	dndl[0] = 0.0;
1137
1138	// Start of outer iteration loop
1139
1140	G4int outerCounter = 0, innerCounter = 0; // three times.
1141	G4bool eliminateThisParticle = true;
1142	G4bool resetEnergies = true;
1143	while( ++outerCounter < 3 )
1144	{
1145	for( l=1; l<20; l++ )
1146	{
1147	xval = (binl[l]+binl[l-1])/2.; // x = lambda /GeV
1148	if( xval > 1./pt )
1149	dndl[l] = dndl[l-1];
1150	else
1151	dndl[l] = dndl[l-1] +
1152	aspar/std::sqrt( std::pow((1.+asparxvalasparxval),3) )
1153	(binl[l]-binl[l-1]) * et /
1154	std::sqrt( ptxvaletptxvalet + ptpt + pvMass*pvMass );
1155	}
1156
1157	// Start of inner iteration loop
1158
1159	innerCounter = 0; // try this not more than 7 times.
1160	while( ++innerCounter < 7 )
1161	{
1162	l = 1;
1163	ran = G4UniformRand()*dndl[19];
1164	while( ( ran >= dndl[l] ) && ( l < 20 ) )l++;
1165	l = Imin( 19, l );
1166	xval = Amin( 1.0, pt(binl[l-1] + G4UniformRand()(binl[l]-binl[l-1]) ) );
1167	if( pv[i].getSide() < 0 ) xval *= -1.;
1168	pv[i].setMomentumAndUpdate( xval*et ); // Set the z-momentum
1169	pvEnergy = pv[i].getEnergy();
1170	if( pv[i].getSide() > 0 ) // Forward side
1171	{
1172	if ( i < 2 )
1173	{
1174	ekin = tavai1 - ekin1;
1175	if (ekin < 0.) ekin = 0.04*std::fabs(normal());
1176	G4double pp1 = pv[i].Length();
1177	if (pp1 >= 1.e-6)
1178	{
1179	G4double pp = std::sqrt(ekin(ekin+2pvMass));
1180	pp = Amax(0., pppp - ptpt);
1181	pp = std::sqrt(pp)*pv[i].getSide()/std::fabs(G4double(pv[i].getSide())); // cast for aCC
1182	pv[i].setMomentumAndUpdate( pp );
1183	}
1184	else
1185	{
1186	pv[i].setMomentum(0.,0.,0.);
1187	pv[i].setKineticEnergyAndUpdate( ekin);
1188	}
1189	pvmx[4].Add( pvmx[4], pv[i]);
1190	outerCounter = 2;
1191	resetEnergies = false;
1192	eliminateThisParticle = false;
1193	break;
1194	}
1195	else if( (ekin1+pvEnergy-pvMass) < 0.95*tavai1 )
1196	{
1197	pvmx[4].Add( pvmx[4], pv[i] );
1198	ekin1 += pvEnergy - pvMass;
1199	pvmx[6].Add( pvmx[4], pvmx[5] );
1200	pvmx[6].setMomentum( 0.0 );
1201	outerCounter = 2; // leave outer loop
1202	eliminateThisParticle = false; // don't eliminate this particle
1203	resetEnergies = false;
1204	break; // next particle
1205	}
1206	if( innerCounter > 5 ) break; // leave inner loop
1207
1208	if( tavai2 >= pvMass )
1209	{ // switch sides
1210	pv[i].setSide( -1 );
1211	tavai1 += pvMass;
1212	tavai2 -= pvMass;
1213	iavai2++;
1214	}
1215	}
1216	else
1217	{ // backward side
1218	xval = Amin(0.999,0.95+0.05*targ/20.0);
1219	if( (ekin2+pvEnergy-pvMass) < xval*tavai2 )
1220	{
1221	pvmx[5].Add( pvmx[5], pv[i] );
1222	ekin2 += pvEnergy - pvMass;
1223	pvmx[6].Add( pvmx[4], pvmx[5] );
1224	pvmx[6].setMomentum( 0.0 ); // set z-momentum
1225	outerCounter = 2; // leave outer iteration
1226	eliminateThisParticle = false; // don't eliminate this particle
1227	resetEnergies = false;
1228	break; // leave inner iteration
1229	}
1230	if( innerCounter > 5 )break; // leave inner iteration
1231
1232	if( tavai1 >= pvMass )
1233	{ // switch sides
1234	pv[i].setSide( 1 );
1235	tavai1 -= pvMass;
1236	tavai2 += pvMass;
1237	iavai2--;
1238	}
1239	}
1240	pv[i].setMomentum( pv[i].getMomentum().x() * 0.9,
1241	pv[i].getMomentum().y() * 0.9);
1242	pt *= 0.9;
1243	dndl[19] *= 0.9;
1244	} // closes inner loop
1245
1246	if (resetEnergies)
1247	{
1248	if (verboseLevel > 1) {
1249	G4cout << " Reset energies for index " << i << " "
1250	<< ekin1 << " " << tavai1 << G4endl;
1251	pv[i].Print(i);
1252	}
1253	ekin1 = 0.0;
1254	ekin2 = 0.0;
1255	pvmx[4].setZero();
1256	pvmx[5].setZero();
1257
1258	for( l=i+1; l < vecLen; l++ )
1259	{
1260	if( (pv[l].getMass() < protonMass) \|\| (pv[l].getSide() > 0) )
1261	{
1262	pvEnergy = pv[l].getMass() + 0.95*pv[l].getKineticEnergy();
1263	pv[l].setEnergyAndUpdate( pvEnergy );
1264	if( pv[l].getSide() > 0)
1265	{
1266	ekin1 += pv[l].getKineticEnergy();
1267	pvmx[4].Add( pvmx[4], pv[l] );
1268	}
1269	else
1270	{
1271	ekin2 += pv[l].getKineticEnergy();
1272	pvmx[5].Add( pvmx[5], pv[l] );
1273	}
1274	}
1275	}
1276	}
1277	} // closes outer iteration
1278
1279	if( eliminateThisParticle ) // not enough energy,
1280	{ // eliminate this particle
1281	if (verboseLevel > 1) {
1282	G4cout << " Eliminate particle index " << i << G4endl;
1283	pv[i].Print(i);
1284	}
1285	if(i != vecLen-1)
1286	{
1287	for( j=i; j < vecLen-1; j++ )
1288	{ // shift down
1289	pv[j] = pv[j+1];
1290	}
1291	}
1292	vecLen--;
1293	if(vecLen < 2) return;
1294	i++;
1295	pvmx[6].Add( pvmx[4], pvmx[5] );
1296	pvmx[6].setMomentum( 0.0 ); // set z-momentum
1297	}
1298	} // closes main for loop
1299	if (verboseLevel > 1)
1300	{ G4cout << " pv Vector after lambda fragmentation " << vecLen << G4endl;
1301	pvI.Print(-1);
1302	pvT.Print(-1);
1303	for (i=0; i < vecLen ; i++) pv[i].Print(i);
1304	for (i=0; i < 10; i++) pvmx[i].Print(i);
1305	}
1306
1307
1308	// Backward nucleons produced with a cluster model
1309
1310	G4double gpar[] = {2.6, 2.6, 1.80, 1.30, 1.20};
1311	G4double cpar[] = {0.6, 0.6, 0.35, 0.15, 0.10};
1312
1313	if (npg > 0)
1314	{
1315	G4double rmg = rmg0;
1316	if (npg > 1)
1317	{
1318	G4int npg1 = npg-1;
1319	if (npg1 >4) npg1 = 4;
1320	rmg += std::pow( -std::log(1.-G4UniformRand()), cpar[npg1])/gpar[npg1];
1321	}
1322	G4double ga = 1.2;
1323	G4double ekit1 = 0.04, ekit2 = 0.6;
1324	if(incidentKineticEnergy < 5.)
1325	{
1326	ekit1 *= sqr(incidentKineticEnergy)/25.;
1327	ekit2 *= sqr(incidentKineticEnergy)/25.;
1328	}
1329	G4double avalue = (1.-ga)/(std::pow(ekit2,1.-ga)-std::pow(ekit1,1.-ga));
1330	for (i = 0; i < vecLen; i++)
1331	{
1332	if (pv[i].getSide() == -3)
1333	{
1334	G4double ekit = std::pow(G4UniformRand()*(1.-ga)/avalue + std::pow(ekit1,1.-ga), 1./(1.-ga) );
1335	G4double cost = Amax(-1., Amin(1., std::log(2.23*G4UniformRand()+0.383)/0.96));
1336	G4double sint = std::sqrt(1. - cost*cost);
1337	G4double phi = twopi*G4UniformRand();
1338	G4double pp = std::sqrt(ekit(ekit+2pv[i].getMass()));
1339	pv[i].setMomentum( ppsintstd::sin(phi),
1340	ppsintstd::cos(phi),
1341	pp*cost );
1342	pv[i].Lor( pv[i], pvmx[2] );
1343	pvmx[5].Add( pvmx[5], pv[i] );
1344	}
1345	}
1346	}
1347
1348	if (vecLen <= 2) {
1349	successful = false;
1350	return;
1351	}
1352
1353	// Lorentz transformation in lab system
1354
1355	targ = 0;
1356	for( i=0; i < vecLen; i++ )
1357	{
1358	if( pv[i].getType() == baryonType )targ++;
1359	if( pv[i].getType() == antiBaryonType )targ--;
1360	if(verboseLevel > 1) pv[i].Print(i);
1361	pv[i].Lor( pv[i], pvmx[1] );
1362	if(verboseLevel > 1) pv[i].Print(i);
1363	}
1364	if ( targ <1) targ = 1;
1365
1366	G4bool dum=0;
1367	if( lead )
1368	{
1369	for( i=0; i<vecLen; i++ )
1370	{
1371	if( pv[i].getCode() == lead )
1372	{
1373	dum = false;
1374	break;
1375	}
1376	}
1377	if( dum )
1378	{
1379	i = 0;
1380
1381	if( ( (leadParticle.getType() == baryonType \|\|
1382	leadParticle.getType() == antiBaryonType)
1383	&& (pv[1].getType() == baryonType \|\|
1384	pv[1].getType() == antiBaryonType))
1385	\|\| ( (leadParticle.getType() == mesonType)
1386	&& (pv[1].getType() == mesonType)))
1387	{
1388	i = 1;
1389	}
1390	ekin = pv[i].getKineticEnergy();
1391	pv[i] = leadParticle;
1392	if( pv[i].getFlag() )
1393	pv[i].setTOF( -1.0 );
1394	else
1395	pv[i].setTOF( 1.0 );
1396	pv[i].setKineticEnergyAndUpdate( ekin );
1397	}
1398	}
1399
1400	pvmx[3].setMass( incidentMass);
1401	pvmx[3].setMomentumAndUpdate( 0.0, 0.0, incidentTotalMomentum );
1402
1403	G4double ekin0 = pvmx[3].getKineticEnergy();
1404
1405	pvmx[4].setMass( protonMass * targ);
1406	pvmx[4].setEnergy( protonMass * targ);
1407	pvmx[4].setKineticEnergy(0.);
1408	pvmx[4].setMomentum(0., 0., 0.);
1409	ekin = pvmx[3].getEnergy() + pvmx[4].getEnergy();
1410
1411	pvmx[5].Add( pvmx[3], pvmx[4] );
1412	pvmx[3].Lor( pvmx[3], pvmx[5] );
1413	pvmx[4].Lor( pvmx[4], pvmx[5] );
1414
1415	G4double tecm = pvmx[3].getEnergy() + pvmx[4].getEnergy();
1416
1417	pvmx[7].setZero();
1418
1419	ekin1 = 0.0;
1420	G4double teta, wgt;
1421
1422	for( i=0; i < vecLen; i++ )
1423	{
1424	pvmx[7].Add( pvmx[7], pv[i] );
1425	ekin1 += pv[i].getKineticEnergy();
1426	ekin -= pv[i].getMass();
1427	}
1428
1429	if( vecLen > 1 && vecLen < 19 )
1430	{
1431	G4bool constantCrossSection = true;
1432	G4HEVector pw[19];
1433	for(i=0; i<vecLen; i++) pw[i] = pv[i];
1434	wgt = NBodyPhaseSpace( tecm, constantCrossSection, pw, vecLen );
1435	ekin = 0.0;
1436	for( i=0; i < vecLen; i++ )
1437	{
1438	pvmx[6].setMass( pw[i].getMass());
1439	pvmx[6].setMomentum( pw[i].getMomentum() );
1440	pvmx[6].SmulAndUpdate( pvmx[6], 1. );
1441	pvmx[6].Lor( pvmx[6], pvmx[4] );
1442	ekin += pvmx[6].getKineticEnergy();
1443	}
1444	teta = pvmx[7].Ang( pvmx[3] );
1445	if (verboseLevel > 1)
1446	G4cout << " vecLen > 1 && vecLen < 19 " << teta << " " << ekin0
1447	<< " " << ekin1 << " " << ekin << G4endl;
1448	}
1449
1450	if( ekin1 != 0.0 )
1451	{
1452	pvmx[6].setZero();
1453	wgt = ekin/ekin1;
1454	ekin1 = 0.;
1455	for( i=0; i < vecLen; i++ )
1456	{
1457	pvMass = pv[i].getMass();
1458	ekin = pv[i].getKineticEnergy() * wgt;
1459	pv[i].setKineticEnergyAndUpdate( ekin );
1460	ekin1 += ekin;
1461	pvmx[6].Add( pvmx[6], pv[i] );
1462	}
1463	teta = pvmx[6].Ang( pvmx[3] );
1464	if (verboseLevel > 1) {
1465	G4cout << " ekin1 != 0 " << teta << " " << ekin0 << " "
1466	<< ekin1 << G4endl;
1467	incidentParticle.Print(0);
1468	targetParticle.Print(1);
1469	for(i=0;i<vecLen;i++) pv[i].Print(i);
1470	}
1471	}
1472
1473	// Do some smearing in the transverse direction due to Fermi motion
1474
1475	G4double ry = G4UniformRand();
1476	G4double rz = G4UniformRand();
1477	G4double rx = twopi*rz;
1478	G4double a1 = std::sqrt(-2.0*std::log(ry));
1479	G4double rantarg1 = a1std::cos(rx)0.02*targ/G4double(vecLen);
1480	G4double rantarg2 = a1std::sin(rx)0.02*targ/G4double(vecLen);
1481
1482	for (i = 0; i < vecLen; i++)
1483	pv[i].setMomentum( pv[i].getMomentum().x()+rantarg1,
1484	pv[i].getMomentum().y()+rantarg2 );
1485
1486	if (verboseLevel > 1) {
1487	pvmx[6].setZero();
1488	for (i = 0; i < vecLen; i++) pvmx[6].Add( pvmx[6], pv[i] );
1489	teta = pvmx[6].Ang( pvmx[3] );
1490	G4cout << " After smearing " << teta << G4endl;
1491	}
1492
1493	// Rotate in the direction of the primary particle momentum (z-axis).
1494	// This does disturb our inclusive distributions somewhat, but it is
1495	// necessary for momentum conservation
1496
1497	// Also subtract binding energies and make some further corrections
1498	// if required
1499
1500	G4double dekin = 0.0;
1501	G4int npions = 0;
1502	G4double ek1 = 0.0;
1503	G4double alekw, xxh;
1504	G4double cfa = 0.025((atomicWeight-1.)/120.)std::exp(-(atomicWeight-1.)/120.);
1505	G4double alem[] = {1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00, 10.00};
1506	G4double val0[] = {0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65, 0.70};
1507
1508	if (verboseLevel > 1)
1509	G4cout << " Rotation in Direction of primary particle (Defs1)" << G4endl;
1510
1511	for (i = 0; i < vecLen; i++)
1512	{
1513	if(verboseLevel > 1) pv[i].Print(i);
1514	pv[i].Defs1( pv[i], pvI );
1515	if(verboseLevel > 1) pv[i].Print(i);
1516	if (atomicWeight > 1.5)
1517	{
1518	ekin = Amax( 1.e-6,pv[i].getKineticEnergy() - cfa( 1. + 0.5normal()));
1519	alekw = std::log( incidentKineticEnergy );
1520	xxh = 1.;
1521	if(incidentCode == pionPlusCode \|\| incidentCode == pionMinusCode)
1522	{
1523	if(pv[i].getCode() == pionZeroCode)
1524	{
1525	if(G4UniformRand() < std::log(atomicWeight))
1526	{
1527	if (alekw > alem[0])
1528	{
1529	for (j = 1; j < 8; j++)
1530	{
1531	if(alekw < alem[j]) break;
1532	}
1533	xxh = (val0[j]-val0[j-1])/(alem[j]-alem[j-1])*alekw
1534	+ val0[j-1] - (val0[j]-val0[j-1])/(alem[j]-alem[j-1])*alem[j-1];
1535	xxh = 1. - xxh;
1536	}
1537	}
1538	}
1539	}
1540	dekin += ekin*(1.-xxh);
1541	ekin *= xxh;
1542	pv[i].setKineticEnergyAndUpdate( ekin );
1543	pvCode = pv[i].getCode();
1544	if ((pvCode == pionPlusCode) \|\| (pvCode == pionMinusCode) \|\| (pvCode == pionZeroCode))
1545	{
1546	npions += 1;
1547	ek1 += ekin;
1548	}
1549	}
1550	}
1551	if( (ek1 > 0.0) && (npions > 0) )
1552	{
1553	dekin = 1.+dekin/ek1;
1554	for (i = 0; i < vecLen; i++)
1555	{
1556	pvCode = pv[i].getCode();
1557	if((pvCode == pionPlusCode) \|\| (pvCode == pionMinusCode) \|\| (pvCode == pionZeroCode))
1558	{
1559	ekin = Amax( 1.0e-6, pv[i].getKineticEnergy() * dekin );
1560	pv[i].setKineticEnergyAndUpdate( ekin );
1561	}
1562	}
1563	}
1564	if (verboseLevel > 1)
1565	{ G4cout << " Lab-System " << ek1 << " " << npions << G4endl;
1566	incidentParticle.Print(0);
1567	targetParticle.Print(1);
1568	for (i=0; i<vecLen; i++) pv[i].Print(i);
1569	}
1570
1571	// Add black track particles
1572	// the total number of particles produced is restricted to 198
1573	// this may have influence on very high energies
1574
1575	if (verboseLevel > 1)
1576	G4cout << " Evaporation : " << atomicWeight << " "
1577	<< excitationEnergyGNP << " " << excitationEnergyDTA << G4endl;
1578
1579	G4double sprob = 0.;
1580	if (incidentKineticEnergy > 5.)
1581	// sprob = Amin(1., (0.394-0.063std::log(atomicWeight))std::log(incidentKineticEnergy-4.) );
1582	sprob = Amin(1., 0.000314atomicWeightstd::log(incidentKineticEnergy-4.));
1583	if( atomicWeight > 1.5 && G4UniformRand() > sprob )
1584	{
1585
1586	G4double cost, sint, pp, eka;
1587	G4int spall(0), nbl(0);
1588
1589	// first add protons and neutrons
1590
1591	if( excitationEnergyGNP >= 0.001 )
1592	{
1593	// nbl = number of proton/neutron black track particles
1594	// tex is their total kinetic energy (GeV)
1595
1596	nbl = Poisson( (1.5+1.25targ)excitationEnergyGNP/
1597	(excitationEnergyGNP+excitationEnergyDTA));
1598	if( targ+nbl > atomicWeight ) nbl = (int)(atomicWeight - targ);
1599	if (verboseLevel > 1)
1600	G4cout << " evaporation " << targ << " " << nbl << " "
1601	<< sprob << G4endl;
1602	spall = targ;
1603	if( nbl > 0)
1604	{
1605	ekin = (excitationEnergyGNP)/nbl;
1606	ekin2 = 0.0;
1607	for( i=0; i<nbl; i++ )
1608	{
1609	if( G4UniformRand() < sprob )
1610	{
1611	if(verboseLevel > 1) G4cout << " Particle skipped " << G4endl;
1612	continue;
1613	}
1614	if( ekin2 > excitationEnergyGNP) break;
1615	ran = G4UniformRand();
1616	ekin1 = -ekinstd::log(ran) - cfa(1.0+0.5*normal());
1617	if (ekin1 < 0) ekin1 = -0.010*std::log(ran);
1618	ekin2 += ekin1;
1619	if( ekin2 > excitationEnergyGNP)
1620	ekin1 = Amax( 1.0e-6, excitationEnergyGNP-(ekin2-ekin1) );
1621	if( G4UniformRand() > (1.0-atomicNumber/(atomicWeight)))
1622	pv[vecLen] = Proton;
1623	else
1624	pv[vecLen] = Neutron;
1625	spall++;
1626	cost = G4UniformRand() * 2.0 - 1.0;
1627	sint = std::sqrt(std::fabs(1.0-cost*cost));
1628	phi = twopi * G4UniformRand();
1629	pv[vecLen].setFlag( true ); // true is the same as IPA(i)<0
1630	pv[vecLen].setSide( -4 );
1631	pv[vecLen].setTOF( 1.0 );
1632	pvMass = pv[vecLen].getMass();
1633	pvEnergy = ekin1 + pvMass;
1634	pp = std::sqrt( std::fabs( pvEnergypvEnergy - pvMasspvMass ) );
1635	pv[vecLen].setMomentumAndUpdate( ppsintstd::sin(phi),
1636	ppsintstd::cos(phi),
1637	pp*cost );
1638	if (verboseLevel > 1) pv[vecLen].Print(vecLen);
1639	vecLen++;
1640	}
1641	if( (atomicWeight >= 10.0 ) && (incidentKineticEnergy <= 2.0) )
1642	{
1643	G4int ika, kk = 0;
1644	eka = incidentKineticEnergy;
1645	if( eka > 1.0 )eka *= eka;
1646	eka = Amax( 0.1, eka );
1647	ika = G4int(3.6std::exp((atomicNumberatomicNumber
1648	/atomicWeight-35.56)/6.45)/eka);
1649	if( ika > 0 )
1650	{
1651	for( i=(vecLen-1); i>=0; i-- )
1652	{
1653	if( (pv[i].getCode() == protonCode) && pv[i].getFlag() )
1654	{
1655	pTemp = pv[i];
1656	pv[i].setDefinition("Neutron");
1657	pv[i].setMomentumAndUpdate(pTemp.getMomentum());
1658	if (verboseLevel > 1) pv[i].Print(i);
1659	if( ++kk > ika ) break;
1660	}
1661	}
1662	}
1663	}
1664	}
1665	}
1666
1667	// finished adding proton/neutron black track particles
1668	// now, try to add deuterons, tritons and alphas
1669
1670	if( excitationEnergyDTA >= 0.001 )
1671	{
1672	nbl = Poisson( (1.5+1.25targ)excitationEnergyDTA
1673	/(excitationEnergyGNP+excitationEnergyDTA));
1674
1675	// nbl is the number of deutrons, tritons, and alphas produced
1676
1677	if (verboseLevel > 1)
1678	G4cout << " evaporation " << targ << " " << nbl << " "
1679	<< sprob << G4endl;
1680	if( nbl > 0 )
1681	{
1682	ekin = excitationEnergyDTA/nbl;
1683	ekin2 = 0.0;
1684	for( i=0; i<nbl; i++ )
1685	{
1686	if( G4UniformRand() < sprob )
1687	{
1688	if(verboseLevel > 1) G4cout << " Particle skipped " << G4endl;
1689	continue;
1690	}
1691	if( ekin2 > excitationEnergyDTA) break;
1692	ran = G4UniformRand();
1693	ekin1 = -ekinstd::log(ran)-cfa(1.+0.5*normal());
1694	if( ekin1 < 0.0 ) ekin1 = -0.010*std::log(ran);
1695	ekin2 += ekin1;
1696	if( ekin2 > excitationEnergyDTA)
1697	ekin1 = Amax( 1.0e-6, excitationEnergyDTA-(ekin2-ekin1));
1698	cost = G4UniformRand()*2.0 - 1.0;
1699	sint = std::sqrt(std::fabs(1.0-cost*cost));
1700	phi = twopi*G4UniformRand();
1701	ran = G4UniformRand();
1702	if( ran <= 0.60 )
1703	pv[vecLen] = Deuteron;
1704	else if (ran <= 0.90)
1705	pv[vecLen] = Triton;
1706	else
1707	pv[vecLen] = Alpha;
1708	spall += (int)(pv[vecLen].getMass() * 1.066);
1709	if( spall > atomicWeight ) break;
1710	pv[vecLen].setFlag( true ); // true is the same as IPA(i)<0
1711	pv[vecLen].setSide( -4 );
1712	pvMass = pv[vecLen].getMass();
1713	pv[vecLen].setTOF( 1.0 );
1714	pvEnergy = pvMass + ekin1;
1715	pp = std::sqrt( std::fabs( pvEnergypvEnergy - pvMasspvMass ) );
1716	pv[vecLen].setMomentumAndUpdate( ppsintstd::sin(phi),
1717	ppsintstd::cos(phi),
1718	pp*cost );
1719	if (verboseLevel > 1) pv[vecLen].Print(vecLen);
1720	vecLen++;
1721	}
1722	}
1723	}
1724	}
1725	if( centerOfMassEnergy <= (4.0+G4UniformRand()) )
1726	{
1727	for( i=0; i<vecLen; i++ )
1728	{
1729	G4double etb = pv[i].getKineticEnergy();
1730	if( etb >= incidentKineticEnergy )
1731	pv[i].setKineticEnergyAndUpdate( incidentKineticEnergy );
1732	}
1733	}
1734
1735	if(verboseLevel > 1)
1736	{ G4cout << "Call TuningOfHighEnergyCacsading vecLen = " << vecLen << G4endl;
1737	incidentParticle.Print(0);
1738	targetParticle.Print(1);
1739	for (i=0; i<vecLen; i++) pv[i].Print(i);
1740	}
1741
1742	TuningOfHighEnergyCascading( pv, vecLen,
1743	incidentParticle, targetParticle,
1744	atomicWeight, atomicNumber);
1745
1746	if(verboseLevel > 1)
1747	{ G4cout << " After Tuning: " << G4endl;
1748	for(i=0; i<vecLen; i++) pv[i].Print(i);
1749	}
1750
1751	// Calculate time delay for nuclear reactions
1752
1753	G4double tof = incidentTOF;
1754	if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0)
1755	&& (incidentKineticEnergy <= 0.2) )
1756	tof -= 500.0 * std::exp(-incidentKineticEnergy /0.04) * std::log( G4UniformRand() );
1757
1758	for(i=0; i<vecLen; i++)
1759	{
1760	if(pv[i].getName() == "KaonZero" \|\| pv[i].getName() == "AntiKaonZero")
1761	{
1762	pvmx[0] = pv[i];
1763	if(G4UniformRand() < 0.5) pv[i].setDefinition("KaonZeroShort");
1764	else pv[i].setDefinition("KaonZeroLong");
1765	pv[i].setMomentumAndUpdate(pvmx[0].getMomentum());
1766	}
1767	}
1768
1769	successful = true;
1770	delete [] pvmx;
1771	G4int testCurr=0;
1772	G4double totKin=0;
1773	for(testCurr=0; testCurr<vecLen; testCurr++)
1774	{
1775	totKin+=pv[testCurr].getKineticEnergy();
1776	if(totKin>incidentKineticEnergy*1.05)
1777	{
1778	vecLen = testCurr;
1779	break;
1780	}
1781	}
1782
1783	return;
1784	}
1785
1786	void
1787	G4HEInelastic::TuningOfHighEnergyCascading(G4HEVector pv[],
1788	G4int &vecLen,
1789	G4HEVector incidentParticle,
1790	G4HEVector targetParticle,
1791	G4double atomicWeight,
1792	G4double atomicNumber)
1793	{
1794	G4int i,j;
1795	G4double incidentKineticEnergy = incidentParticle.getKineticEnergy();
1796	G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
1797	G4double incidentCharge = incidentParticle.getCharge();
1798	G4double incidentMass = incidentParticle.getMass();
1799	G4double targetMass = targetParticle.getMass();
1800	G4int pionPlusCode = PionPlus.getCode();
1801	G4int pionMinusCode = PionMinus.getCode();
1802	G4int pionZeroCode = PionZero.getCode();
1803	G4int protonCode = Proton.getCode();
1804	G4int neutronCode = Neutron.getCode();
1805	G4HEVector *pvmx = new G4HEVector [10];
1806	G4double *reddec = new G4double [7];
1807
1808	if (incidentKineticEnergy > (25.+G4UniformRand()*75.) ) {
1809	G4double reden = -0.7 + 0.29*std::log10(incidentKineticEnergy);
1810	// G4double redat = 1.0 - 0.40*std::log10(atomicWeight);
1811	// G4double redat = 0.5 - 0.18*std::log10(atomicWeight);
1812	G4double redat = 0.722 - 0.278*std::log10(atomicWeight);
1813	G4double pmax = -200.;
1814	G4double pmapim = -200.;
1815	G4double pmapi0 = -200.;
1816	G4double pmapip = -200.;
1817	G4int ipmax = 0;
1818	G4int maxpim = 0;
1819	G4int maxpi0 = 0;
1820	G4int maxpip = 0;
1821	G4int iphmf;
1822	if ( (G4UniformRand() > (atomicWeight/100.-0.28))
1823	&& (std::fabs(incidentCharge) > 0.) )
1824	{
1825	for (i=0; i < vecLen; i++)
1826	{
1827	iphmf = pv[i].getCode();
1828	G4double ppp = pv[i].Length();
1829	if ( ppp > pmax)
1830	{
1831	pmax = ppp; ipmax = i;
1832	}
1833	if (iphmf == pionPlusCode && ppp > pmapip )
1834	{
1835	pmapip = ppp; maxpip = i;
1836	}
1837	else if (iphmf == pionZeroCode && ppp > pmapi0)
1838	{
1839	pmapi0 = ppp; maxpi0 = i;
1840	}
1841	else if (iphmf == pionMinusCode && ppp > pmapim)
1842	{
1843	pmapim = ppp; maxpim = i;
1844	}
1845	}
1846	}
1847	if(verboseLevel > 1)
1848	{
1849	G4cout << "ipmax, pmax " << ipmax << " " << pmax << G4endl;
1850	G4cout << "maxpip,pmapip " << maxpip << " " << pmapip << G4endl;
1851	G4cout << "maxpi0,pmapi0 " << maxpi0 << " " << pmapi0 << G4endl;
1852	G4cout << "maxpim,pmapim " << maxpim << " " << pmapim << G4endl;
1853	}
1854
1855	if ( vecLen > 2)
1856	{
1857	for (i=2; i<vecLen; i++)
1858	{
1859	iphmf = pv[i].getCode();
1860	if ( ((iphmf==protonCode)\|\|(iphmf==neutronCode)\|\|(pv[i].getType()=="Nucleus"))
1861	&& (pv[i].Length()<1.5)
1862	&& ((G4UniformRand()<reden)\|\|(G4UniformRand()<redat)))
1863	{
1864	pv[i].setMomentumAndUpdate( 0., 0., 0.);
1865	if(verboseLevel > 1)
1866	{
1867	G4cout << "zero Momentum for particle " << G4endl;
1868	pv[i].Print(i);
1869	}
1870	}
1871	}
1872	}
1873	if (maxpi0 == ipmax)
1874	{
1875	if (G4UniformRand() < pmapi0/incidentTotalMomentum)
1876	{
1877	if ((incidentCharge > 0.5) && (maxpip != 0))
1878	{
1879	G4ParticleMomentum mompi0 = pv[maxpi0].getMomentum();
1880	pv[maxpi0].setMomentumAndUpdate( pv[maxpip].getMomentum() );
1881	pv[maxpip].setMomentumAndUpdate( mompi0);
1882	if(verboseLevel > 1)
1883	{
1884	G4cout << " exchange Momentum for " << maxpi0 << " and " << maxpip << G4endl;
1885	}
1886	}
1887	else if ((incidentCharge < -0.5) && (maxpim != 0))
1888	{
1889	G4ParticleMomentum mompi0 = pv[maxpi0].getMomentum();
1890	pv[maxpi0].setMomentumAndUpdate( pv[maxpim].getMomentum() );
1891	pv[maxpim].setMomentumAndUpdate( mompi0);
1892	if(verboseLevel > 1)
1893	{
1894	G4cout << " exchange Momentum for " << maxpi0 << " and " << maxpip << G4endl;
1895	}
1896	}
1897	}
1898	}
1899	G4double bntot = - incidentParticle.getBaryonNumber() - atomicWeight;
1900	for (i=0; i<vecLen; i++) bntot += pv[i].getBaryonNumber();
1901	if(atomicWeight < 1.5) { bntot = 0.; }
1902	else { bntot = 1. + bntot/atomicWeight; }
1903	if(atomicWeight > (75.+G4UniformRand()*50.)) bntot = 0.;
1904	if(verboseLevel > 1)
1905	{
1906	G4cout << " Calculated Baryon- Number " << bntot << G4endl;
1907	}
1908
1909	j = 0;
1910	for (i=0; i<vecLen; i++)
1911	{
1912	G4double ppp = pv[i].Length();
1913	if ( ppp > 1.e-6)
1914	{
1915	iphmf = pv[i].getCode();
1916	if( (bntot > 0.3)
1917	&& ((iphmf == protonCode) \|\| (iphmf == neutronCode)
1918	\|\| (pv[i].getType() == "Nucleus") )
1919	&& (G4UniformRand() < 0.25)
1920	&& (ppp < 1.2) )
1921	{
1922	if(verboseLevel > 1)
1923	{
1924	G4cout << " skip Baryon: " << G4endl;
1925	pv[i].Print(i);
1926	}
1927	continue;
1928
1929	}
1930	if (j != i)
1931	{
1932	pv[j] = pv[i];
1933	}
1934	j++;
1935	}
1936	}
1937	vecLen = j;
1938	}
1939
1940	pvmx[0] = incidentParticle;
1941	pvmx[1] = targetParticle;
1942	pvmx[8].setZero();
1943	pvmx[2].Add(pvmx[0], pvmx[1]);
1944	pvmx[3].Lor(pvmx[0], pvmx[2]);
1945	pvmx[4].Lor(pvmx[1], pvmx[2]);
1946
1947	if (verboseLevel > 1) {
1948	pvmx[0].Print(0);
1949	incidentParticle.Print(0);
1950	pvmx[1].Print(1);
1951	targetParticle.Print(1);
1952	pvmx[2].Print(2);
1953	pvmx[3].Print(3);
1954	pvmx[4].Print(4);
1955	}
1956
1957	// Calculate leading particle effect in which a single final state
1958	// particle carries away nearly the maximum allowed momentum, while
1959	// all other secondaries have reduced momentum. A secondary is
1960	// proportionately less likely to be a leading particle as the
1961	// difference of its quantum numbers with the primary increases.
1962
1963	G4int ledpar = -1;
1964	G4double redpar = 0.;
1965	G4int incidentS = incidentParticle.getStrangenessNumber();
1966	if (incidentParticle.getName() == "KaonZeroShort" \|\|
1967	incidentParticle.getName() == "KaonZeroLong") {
1968	if(G4UniformRand() < 0.5) {
1969	incidentS = 1;
1970	} else {
1971	incidentS = -1;
1972	}
1973	}
1974
1975	G4int incidentB = incidentParticle.getBaryonNumber();
1976
1977	for (i=0; i<vecLen; i++) {
1978	G4int iphmf = pv[i].getCode();
1979	G4double ppp = pv[i].Length();
1980
1981	if (ppp > 1.e-3) {
1982	pvmx[5].Lor( pv[i], pvmx[2] ); // secondary in CM frame
1983	G4double cost = pvmx[3].CosAng( pvmx[5] );
1984
1985	// For each secondary, find the sum of the differences of its
1986	// quantum numbers with that of the incident particle
1987	// (dM + dQ + dS + dB)
1988
1989	G4int particleS = pv[i].getStrangenessNumber();
1990
1991	if (pv[i].getName() == "KaonZeroShort" \|\|
1992	pv[i].getName() == "KaonZeroLong") {
1993	if (G4UniformRand() < 0.5) {
1994	particleS = 1;
1995	} else {
1996	particleS = -1;
1997	}
1998	}
1999	G4int particleB = pv[i].getBaryonNumber();
2000	G4double hfmass;
2001	if (cost > 0.) {
2002	reddec[0] = std::fabs( incidentMass - pv[i].getMass() );
2003	reddec[1] = std::fabs( incidentCharge - pv[i].getCharge());
2004	reddec[2] = std::fabs( G4double(incidentS - particleS) ); // cast for aCC
2005	reddec[3] = std::fabs( G4double(incidentB - particleB) ); // cast for aCC
2006	hfmass = incidentMass;
2007
2008	} else {
2009	reddec[0] = std::fabs( targetMass - pv[i].getMass() );
2010	reddec[1] = std::fabs( atomicNumber/atomicWeight - pv[i].getCharge());
2011	reddec[2] = std::fabs( G4double(particleS) ); // cast for aCC
2012	reddec[3] = std::fabs( 1. - particleB );
2013	hfmass = targetMass;
2014	}
2015
2016	reddec[5] = reddec[0]+reddec[1]+reddec[2]+reddec[3];
2017	G4double sbqwgt = reddec[5];
2018	if (hfmass < 0.2) {
2019	sbqwgt = (sbqwgt-2.5)*0.10;
2020	if(pv[i].getCode() == pionZeroCode) sbqwgt = 0.15;
2021	} else if (hfmass < 0.6) {
2022	sbqwgt = (sbqwgt-3.0)*0.10;
2023	} else {
2024	sbqwgt = (sbqwgt-2.0)*0.10;
2025	if(pv[i].getCode() == pionZeroCode) sbqwgt = 0.15;
2026	}
2027
2028	ppp = pvmx[5].Length();
2029
2030	// Reduce the longitudinal momentum of the secondary by a factor
2031	// which is a function of the sum of the differences
2032
2033	if (sbqwgt > 0. && ppp > 1.e-6) {
2034	G4double pthmf = pppstd::sqrt(1.-costcost);
2035	G4double plhmf = pppcost(1.-sbqwgt);
2036	pvmx[7].Cross( pvmx[3], pvmx[5] );
2037	pvmx[7].Cross( pvmx[7], pvmx[3] );
2038
2039	if (pvmx[3].Length() > 0.)
2040	pvmx[6].SmulAndUpdate( pvmx[3], plhmf/pvmx[3].Length() );
2041	else if(verboseLevel > 1)
2042	G4cout << "NaNQ in Tuning of High Energy Hadronic Interactions" << G4endl;
2043
2044	if (pvmx[7].Length() > 0.)
2045	pvmx[7].SmulAndUpdate( pvmx[7], pthmf/pvmx[7].Length() );
2046	else if(verboseLevel > 1)
2047	G4cout << "NaNQ in Tuning of High Energy Hadronic Interactions" << G4endl;
2048
2049	pvmx[5].Add3(pvmx[6], pvmx[7] );
2050	pvmx[5].setEnergy( std::sqrt(sqr(pvmx[5].Length()) + sqr(pv[i].getMass())));
2051	pv[i].Lor( pvmx[5], pvmx[4] );
2052	if (verboseLevel > 1) {
2053	G4cout << " Particle Momentum changed to: " << G4endl;
2054	pv[i].Print(i);
2055	}
2056	}
2057
2058	// Choose leading particle
2059	// Neither pi0s, backward nucleons from intra-nuclear cascade,
2060	// nor evaporation fragments can be leading particles
2061
2062	G4int ss = -3;
2063	if (incidentS != 0) ss = 0;
2064	if (iphmf != pionZeroCode && pv[i].getSide() > ss) {
2065	pvmx[7].Sub3( incidentParticle, pv[i] );
2066	reddec[4] = pvmx[7].Length()/incidentTotalMomentum;
2067	reddec[6] = reddec[4]*2./3. + reddec[5]/12.;
2068	reddec[6] = Amax(0., 1. - reddec[6]);
2069	if ( (reddec[5] <= 3.75) && (reddec[6] > redpar) ) {
2070	ledpar = i;
2071	redpar = reddec[6];
2072	}
2073	}
2074	}
2075	pvmx[8].Add3(pvmx[8], pv[i] );
2076	}
2077
2078	if(false) if (ledpar >= 0)
2079	{
2080	if(verboseLevel > 1)
2081	{
2082	G4cout << " Leading Particle found : " << ledpar << G4endl;
2083	pv[ledpar].Print(ledpar);
2084	pvmx[8].Print(-2);
2085	incidentParticle.Print(-1);
2086	}
2087	pvmx[4].Sub3(incidentParticle,pvmx[8]);
2088	pvmx[5].Smul(incidentParticle, incidentParticle.CosAng(pvmx[4])
2089	*pvmx[4].Length()/incidentParticle.Length());
2090	pv[ledpar].Add3(pv[ledpar],pvmx[5]);
2091
2092	pv[ledpar].SmulAndUpdate( pv[ledpar], 1.);
2093	if(verboseLevel > 1)
2094	{
2095	pv[ledpar].Print(ledpar);
2096	}
2097	}
2098
2099	if (conserveEnergy) {
2100	G4double ekinhf = 0.;
2101	for (i=0; i<vecLen; i++) {
2102	ekinhf += pv[i].getKineticEnergy();
2103	if(pv[i].getMass() < 0.7) ekinhf += pv[i].getMass();
2104	}
2105	if(incidentParticle.getMass() < 0.7) ekinhf -= incidentParticle.getMass();
2106
2107	if(ledpar < 0) { // no leading particle chosen
2108	ekinhf = incidentParticle.getKineticEnergy()/ekinhf;
2109	for (i=0; i<vecLen; i++)
2110	pv[i].setKineticEnergyAndUpdate(ekinhf*pv[i].getKineticEnergy());
2111
2112	} else {
2113	// take the energy removed from non-leading particles and
2114	// give it to the leading particle
2115	ekinhf = incidentParticle.getKineticEnergy() - ekinhf;
2116	ekinhf += pv[ledpar].getKineticEnergy();
2117	if(ekinhf < 0.) ekinhf = 0.;
2118	pv[ledpar].setKineticEnergyAndUpdate(ekinhf);
2119	}
2120	}
2121
2122	delete [] reddec;
2123	delete [] pvmx;
2124
2125	return;
2126	}
2127
2128	void
2129	G4HEInelastic::HighEnergyClusterProduction(G4bool &successful,
2130	G4HEVector pv[],
2131	G4int &vecLen,
2132	G4double &excitationEnergyGNP,
2133	G4double &excitationEnergyDTA,
2134	G4HEVector incidentParticle,
2135	G4HEVector targetParticle,
2136	G4double atomicWeight,
2137	G4double atomicNumber)
2138	{
2139	// For low multiplicity in the first intranuclear interaction the cascading process
2140	// as described in G4HEInelastic::MediumEnergyCascading does not work
2141	// satisfactorily. From experimental data it is strongly suggested to use
2142	// a two- body resonance model.
2143	//
2144	// All quantities on the G4HEVector Array pv are in GeV- units.
2145
2146	G4int protonCode = Proton.getCode();
2147	G4double protonMass = Proton.getMass();
2148	G4int neutronCode = Neutron.getCode();
2149	G4double kaonPlusMass = KaonPlus.getMass();
2150	G4int pionPlusCode = PionPlus.getCode();
2151	G4int pionZeroCode = PionZero.getCode();
2152	G4int pionMinusCode = PionMinus.getCode();
2153	G4String mesonType = PionPlus.getType();
2154	G4String baryonType = Proton.getType();
2155	G4String antiBaryonType= AntiProton.getType();
2156
2157	G4double targetMass = targetParticle.getMass();
2158
2159	G4int incidentCode = incidentParticle.getCode();
2160	G4double incidentMass = incidentParticle.getMass();
2161	G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
2162	G4double incidentEnergy = incidentParticle.getEnergy();
2163	G4double incidentKineticEnergy = incidentParticle.getKineticEnergy();
2164	G4String incidentType = incidentParticle.getType();
2165	// G4double incidentTOF = incidentParticle.getTOF();
2166	G4double incidentTOF = 0.;
2167
2168	// some local variables
2169
2170	G4int i, j;
2171
2172	if(verboseLevel > 1) G4cout << " G4HEInelastic::HighEnergyClusterProduction " << G4endl;
2173
2174	successful = false;
2175	if(incidentTotalMomentum < 25. + G4UniformRand()*25.) return;
2176
2177	G4double centerOfMassEnergy = std::sqrt( sqr(incidentMass) + sqr(targetMass)
2178	+2.targetMassincidentEnergy);
2179
2180	G4HEVector pvI = incidentParticle; // for the incident particle
2181	pvI.setSide( 1 );
2182
2183	G4HEVector pvT = targetParticle; // for the target particle
2184	pvT.setMomentumAndUpdate( 0.0, 0.0, 0.0 );
2185	pvT.setSide( -1 );
2186	pvT.setTOF( -1.);
2187	// distribute particles in forward and backward
2188	// hemispheres. Note that only low multiplicity
2189	// events from FirstIntInNuc.... should go into
2190	// this routine.
2191	G4int targ = 0;
2192	G4int ifor = 0;
2193	G4int iback = 0;
2194	G4int pvCode;
2195	G4double pvMass, pvEnergy;
2196
2197	pv[0].setSide( 1 );
2198	pv[1].setSide( -1 );
2199	for(i = 0; i < vecLen; i++)
2200	{
2201	if (i > 1)
2202	{
2203	if( G4UniformRand() < 0.5)
2204	{
2205	pv[i].setSide( 1 );
2206	if (++ifor > 18)
2207	{
2208	pv[i].setSide( -1 );
2209	ifor--;
2210	iback++;
2211	}
2212	}
2213	else
2214	{
2215	pv[i].setSide( -1 );
2216	if (++iback > 18)
2217	{
2218	pv[i].setSide( 1 );
2219	ifor++;
2220	iback--;
2221	}
2222	}
2223	}
2224
2225	pvCode = pv[i].getCode();
2226
2227	if ( ( (incidentCode == protonCode) \|\| (incidentCode == neutronCode)
2228	\|\| (incidentType == mesonType) )
2229	&& ( (pvCode == pionPlusCode) \|\| (pvCode == pionMinusCode) )
2230	&& ( (G4UniformRand() < (10.-incidentTotalMomentum)/6.) )
2231	&& ( (G4UniformRand() < atomicWeight/300.) ) )
2232	{
2233	if (G4UniformRand() > atomicNumber/atomicWeight)
2234	pv[i].setDefinition( "Neutron" );
2235	else
2236	pv[i].setDefinition( "Proton" );
2237	targ++;
2238	}
2239	pv[i].setTOF( incidentTOF );
2240	}
2241	G4double tb = 2. * iback;
2242	if (centerOfMassEnergy < (2+G4UniformRand())) tb = (2.*iback + vecLen)/2.;
2243
2244	G4double nucsup[] = { 1.0, 0.7, 0.5, 0.4, 0.35, 0.3};
2245	G4double psup[] = { 3. , 6. , 20., 50., 100.,1000.};
2246	G4double s = centerOfMassEnergy*centerOfMassEnergy;
2247
2248	G4double xtarg = Amax(0.01, Amin(0.50, 0.312+0.2*std::log(std::log(s))+std::pow(s,1.5)/6000.)
2249	* (std::pow(atomicWeight,0.33)-1.) * tb);
2250	G4int momentumBin = 0;
2251	while( (momentumBin < 6) && (incidentTotalMomentum > psup[momentumBin])) momentumBin++;
2252	momentumBin = Imin(5, momentumBin);
2253	G4double xpnhmf = Amax(0.01,xtarg*nucsup[momentumBin]);
2254	G4double xshhmf = Amax(0.01,xtarg-xpnhmf);
2255	G4double rshhmf = 0.25*xshhmf;
2256	G4double rpnhmf = 0.81*xpnhmf;
2257	G4double xhmf;
2258	G4int nshhmf, npnhmf;
2259	if (rshhmf > 1.1)
2260	{
2261	rshhmf = xshhmf/(rshhmf-1.);
2262	if (rshhmf <= 20.)
2263	xhmf = GammaRand( rshhmf );
2264	else
2265	xhmf = Erlang( G4int(rshhmf+0.5) );
2266	xshhmf *= xhmf/rshhmf;
2267	}
2268	nshhmf = Poisson( xshhmf );
2269	if (rpnhmf > 1.1)
2270	{
2271	rpnhmf = xpnhmf/(rpnhmf -1.);
2272	if (rpnhmf <= 20.)
2273	xhmf = GammaRand( rpnhmf );
2274	else
2275	xhmf = Erlang( G4int(rpnhmf+0.5) );
2276	xpnhmf *= xhmf/rpnhmf;
2277	}
2278	npnhmf = Poisson( xpnhmf );
2279
2280	while (npnhmf > 0)
2281	{
2282	if ( G4UniformRand() > (1. - atomicNumber/atomicWeight))
2283	pv[vecLen].setDefinition( "Proton" );
2284	else
2285	pv[vecLen].setDefinition( "Neutron" );
2286	targ++;
2287	pv[vecLen].setSide( -2 );
2288	pv[vecLen].setFlag( true );
2289	pv[vecLen].setTOF( incidentTOF );
2290	vecLen++;
2291	npnhmf--;
2292	}
2293	while (nshhmf > 0)
2294	{
2295	G4double ran = G4UniformRand();
2296	if (ran < 0.333333 )
2297	pv[vecLen].setDefinition( "PionPlus" );
2298	else if (ran < 0.6667)
2299	pv[vecLen].setDefinition( "PionZero" );
2300	else
2301	pv[vecLen].setDefinition( "PionMinus" );
2302	pv[vecLen].setSide( -2 );
2303	pv[vecLen].setFlag( true );
2304	pv[vecLen].setTOF( incidentTOF );
2305	vecLen++;
2306	nshhmf--;
2307	}
2308
2309	// Mark leading particles for incident strange particles
2310	// and antibaryons, for all other we assume that the first
2311	// and second particle are the leading particles.
2312	// We need this later for kinematic aspects of strangeness conservation.
2313
2314	G4int lead = 0;
2315	G4HEVector leadParticle;
2316	if( (incidentMass >= kaonPlusMass-0.05) && (incidentCode != protonCode)
2317	&& (incidentCode != neutronCode) )
2318	{
2319	G4double pMass = pv[0].getMass();
2320	G4int pCode = pv[0].getCode();
2321	if( (pMass >= kaonPlusMass-0.05) && (pCode != protonCode)
2322	&& (pCode != neutronCode) )
2323	{
2324	lead = pCode;
2325	leadParticle = pv[0];
2326	}
2327	else
2328	{
2329	pMass = pv[1].getMass();
2330	pCode = pv[1].getCode();
2331	if( (pMass >= kaonPlusMass-0.05) && (pCode != protonCode)
2332	&& (pCode != neutronCode) )
2333	{
2334	lead = pCode;
2335	leadParticle = pv[1];
2336	}
2337	}
2338	}
2339
2340	if (verboseLevel > 1)
2341	{ G4cout << " pv Vector after initialization " << vecLen << G4endl;
2342	pvI.Print(-1);
2343	pvT.Print(-1);
2344	for (i=0; i < vecLen ; i++) pv[i].Print(i);
2345	}
2346
2347	G4double tavai = 0.;
2348	for(i=0;i<vecLen;i++) if(pv[i].getSide() != -2) tavai += pv[i].getMass();
2349
2350	while (tavai > centerOfMassEnergy)
2351	{
2352	for (i=vecLen-1; i >= 0; i--)
2353	{
2354	if (pv[i].getSide() != -2)
2355	{
2356	tavai -= pv[i].getMass();
2357	if( i != vecLen-1)
2358	{
2359	for (j=i; j < vecLen; j++)
2360	{
2361	pv[j] = pv[j+1];
2362	}
2363	}
2364	if ( --vecLen < 2)
2365	{
2366	successful = false;
2367	return;
2368	}
2369	break;
2370	}
2371	}
2372	}
2373
2374	// Now produce 3 Clusters:
2375	// 1. forward cluster
2376	// 2. backward meson cluster
2377	// 3. backward nucleon cluster
2378
2379	G4double rmc0 = 0., rmd0 = 0., rme0 = 0.;
2380	G4int ntc = 0, ntd = 0, nte = 0;
2381
2382	for (i=0; i < vecLen; i++)
2383	{
2384	if(pv[i].getSide() > 0)
2385	{
2386	if(ntc < 17)
2387	{
2388	rmc0 += pv[i].getMass();
2389	ntc++;
2390	}
2391	else
2392	{
2393	if(ntd < 17)
2394	{
2395	pv[i].setSide(-1);
2396	rmd0 += pv[i].getMass();
2397	ntd++;
2398	}
2399	else
2400	{
2401	pv[i].setSide(-2);
2402	rme0 += pv[i].getMass();
2403	nte++;
2404	}
2405	}
2406	}
2407	else if (pv[i].getSide() == -1)
2408	{
2409	if(ntd < 17)
2410	{
2411	rmd0 += pv[i].getMass();
2412	ntd++;
2413	}
2414	else
2415	{
2416	pv[i].setSide(-2);
2417	rme0 += pv[i].getMass();
2418	nte++;
2419	}
2420	}
2421	else
2422	{
2423	rme0 += pv[i].getMass();
2424	nte++;
2425	}
2426	}
2427
2428	G4double cpar[] = {0.6, 0.6, 0.35, 0.15, 0.10};
2429	G4double gpar[] = {2.6, 2.6, 1.80, 1.30, 1.20};
2430
2431	G4double rmc = rmc0, rmd = rmd0, rme = rme0;
2432	G4int ntc1 = Imin(4,ntc-1);
2433	G4int ntd1 = Imin(4,ntd-1);
2434	G4int nte1 = Imin(4,nte-1);
2435	if (ntc > 1) rmc = rmc0 + std::pow(-std::log(1.-G4UniformRand()),cpar[ntc1])/gpar[ntc1];
2436	if (ntd > 1) rmd = rmd0 + std::pow(-std::log(1.-G4UniformRand()),cpar[ntd1])/gpar[ntd1];
2437	if (nte > 1) rme = rme0 + std::pow(-std::log(1.-G4UniformRand()),cpar[nte1])/gpar[nte1];
2438	while( (rmc+rmd) > centerOfMassEnergy)
2439	{
2440	if ((rmc == rmc0) && (rmd == rmd0))
2441	{
2442	rmd = 0.999centerOfMassEnergy/(rmc+rmd);
2443	rmc = 0.999centerOfMassEnergy/(rmc+rmd);
2444	}
2445	else
2446	{
2447	rmc = 0.1rmc0 + 0.9rmc;
2448	rmd = 0.1rmd0 + 0.9rmd;
2449	}
2450	}
2451	if(verboseLevel > 1)
2452	G4cout << " Cluster Masses: " << ntc << " " << rmc << " " << ntd
2453	<< " " << rmd << " " << nte << " " << rme << G4endl;
2454
2455	G4HEVector* pvmx = new G4HEVector[11];
2456
2457	pvmx[1].setMass( incidentMass);
2458	pvmx[1].setMomentumAndUpdate( 0., 0., incidentTotalMomentum);
2459	pvmx[2].setMass( targetMass);
2460	pvmx[2].setMomentumAndUpdate( 0., 0., 0.);
2461	pvmx[0].Add( pvmx[1], pvmx[2] );
2462	pvmx[1].Lor( pvmx[1], pvmx[0] );
2463	pvmx[2].Lor( pvmx[2], pvmx[0] );
2464
2465	G4double pf = std::sqrt(Amax(0.0001, sqr(sqr(centerOfMassEnergy) + rmdrmd -rmcrmc)
2466	- 4sqr(centerOfMassEnergy)rmdrmd))/(2.centerOfMassEnergy);
2467	pvmx[3].setMass( rmc );
2468	pvmx[4].setMass( rmd );
2469	pvmx[3].setEnergy( std::sqrt(pfpf + rmcrmc) );
2470	pvmx[4].setEnergy( std::sqrt(pfpf + rmdrmd) );
2471
2472	G4double tvalue = -MAXFLOAT;
2473	G4double bvalue = Amax(0.01, 4.0 + 1.6*std::log(incidentTotalMomentum));
2474	if (bvalue != 0.0) tvalue = std::log(G4UniformRand())/bvalue;
2475	G4double pin = pvmx[1].Length();
2476	G4double tacmin = sqr( pvmx[1].getEnergy() - pvmx[3].getEnergy()) - sqr( pin - pf);
2477	G4double ctet = Amax(-1., Amin(1., 1.+2.(tvalue-tacmin)/Amax(1.e-10, 4.pin*pf)));
2478	G4double stet = std::sqrt(Amax(0., 1.0 - ctet*ctet));
2479	G4double phi = twopi * G4UniformRand();
2480	pvmx[3].setMomentum( pf * stet * std::sin(phi),
2481	pf * stet * std::cos(phi),
2482	pf * ctet );
2483	pvmx[4].Smul( pvmx[3], -1.);
2484
2485	if (nte > 0)
2486	{
2487	G4double ekit1 = 0.04;
2488	G4double ekit2 = 0.6;
2489	G4double gaval = 1.2;
2490	if (incidentKineticEnergy <= 5.)
2491	{
2492	ekit1 *= sqr(incidentKineticEnergy)/25.;
2493	ekit2 *= sqr(incidentKineticEnergy)/25.;
2494	}
2495	G4double avalue = (1.-gaval)/(std::pow(ekit2, 1.-gaval)-std::pow(ekit1, 1.-gaval));
2496	for (i=0; i < vecLen; i++)
2497	{
2498	if (pv[i].getSide() == -2)
2499	{
2500	G4double ekit = std::pow(G4UniformRand()*(1.-gaval)/avalue +std::pow(ekit1, 1.-gaval),
2501	1./(1.-gaval));
2502	pv[i].setKineticEnergyAndUpdate( ekit );
2503	ctet = Amax(-1., Amin(1., std::log(2.23*G4UniformRand()+0.383)/0.96));
2504	stet = std::sqrt( Amax( 0.0, 1. - ctet*ctet ));
2505	phi = G4UniformRand()*twopi;
2506	G4double pp = pv[i].Length();
2507	pv[i].setMomentum( pp * stet * std::sin(phi),
2508	pp * stet * std::cos(phi),
2509	pp * ctet );
2510	pv[i].Lor( pv[i], pvmx[0] );
2511	}
2512	}
2513	}
2514	// pvmx[1] = pvmx[3];
2515	// pvmx[2] = pvmx[4];
2516	pvmx[5].SmulAndUpdate( pvmx[3], -1.);
2517	pvmx[6].SmulAndUpdate( pvmx[4], -1.);
2518
2519	if (verboseLevel > 1) {
2520	G4cout << " General vectors before Phase space Generation " << G4endl;
2521	for (i=0; i<7; i++) pvmx[i].Print(i);
2522	}
2523
2524	G4HEVector* tempV = new G4HEVector[18];
2525	G4bool constantCrossSection = true;
2526	G4double wgt;
2527	G4int npg;
2528
2529	if (ntc > 1)
2530	{
2531	npg = 0;
2532	for (i=0; i < vecLen; i++)
2533	{
2534	if (pv[i].getSide() > 0)
2535	{
2536	tempV[npg++] = pv[i];
2537	}
2538	}
2539	wgt = NBodyPhaseSpace( pvmx[3].getMass(), constantCrossSection, tempV, npg);
2540
2541	npg = 0;
2542	for (i=0; i < vecLen; i++)
2543	{
2544	if (pv[i].getSide() > 0)
2545	{
2546	pv[i].setMomentum( tempV[npg++].getMomentum());
2547	pv[i].SmulAndUpdate( pv[i], 1. );
2548	pv[i].Lor( pv[i], pvmx[5] );
2549	}
2550	}
2551	}
2552	else if(ntc == 1)
2553	{
2554	for(i=0; i<vecLen; i++)
2555	{
2556	if(pv[i].getSide() > 0) pv[i].setMomentumAndUpdate(pvmx[3].getMomentum());
2557	}
2558	}
2559	else
2560	{
2561	}
2562
2563	if (ntd > 1)
2564	{
2565	npg = 0;
2566	for (i=0; i < vecLen; i++)
2567	{
2568	if (pv[i].getSide() == -1)
2569	{
2570	tempV[npg++] = pv[i];
2571	}
2572	}
2573	wgt = NBodyPhaseSpace( pvmx[4].getMass(), constantCrossSection, tempV, npg);
2574
2575	npg = 0;
2576	for (i=0; i < vecLen; i++)
2577	{
2578	if (pv[i].getSide() == -1)
2579	{
2580	pv[i].setMomentum( tempV[npg++].getMomentum());
2581	pv[i].SmulAndUpdate( pv[i], 1.);
2582	pv[i].Lor( pv[i], pvmx[6] );
2583	}
2584	}
2585	}
2586	else if(ntd == 1)
2587	{
2588	for(i=0; i<vecLen; i++)
2589	{
2590	if(pv[i].getSide() == -1) pv[i].setMomentumAndUpdate(pvmx[4].getMomentum());
2591	}
2592	}
2593	else
2594	{
2595	}
2596
2597	if(verboseLevel > 1)
2598	{
2599	G4cout << " Vectors after PhaseSpace generation " << G4endl;
2600	for(i=0; i<vecLen; i++) pv[i].Print(i);
2601	}
2602
2603	// Lorentz transformation in lab system
2604
2605	targ = 0;
2606	for( i=0; i < vecLen; i++ )
2607	{
2608	if( pv[i].getType() == baryonType )targ++;
2609	if( pv[i].getType() == antiBaryonType )targ--;
2610	pv[i].Lor( pv[i], pvmx[2] );
2611	}
2612	if (targ<1) targ = 1;
2613
2614	if(verboseLevel > 1) {
2615	G4cout << " Transformation in Lab- System " << G4endl;
2616	for(i=0; i<vecLen; i++) pv[i].Print(i);
2617	}
2618
2619	G4bool dum(0);
2620	G4double ekin, teta;
2621
2622	if( lead )
2623	{
2624	for( i=0; i<vecLen; i++ )
2625	{
2626	if( pv[i].getCode() == lead )
2627	{
2628	dum = false;
2629	break;
2630	}
2631	}
2632	if( dum )
2633	{
2634	i = 0;
2635
2636	if( ( (leadParticle.getType() == baryonType \|\|
2637	leadParticle.getType() == antiBaryonType)
2638	&& (pv[1].getType() == baryonType \|\|
2639	pv[1].getType() == antiBaryonType))
2640	\|\| ( (leadParticle.getType() == mesonType)
2641	&& (pv[1].getType() == mesonType)))
2642	{
2643	i = 1;
2644	}
2645
2646	ekin = pv[i].getKineticEnergy();
2647	pv[i] = leadParticle;
2648	if( pv[i].getFlag() )
2649	pv[i].setTOF( -1.0 );
2650	else
2651	pv[i].setTOF( 1.0 );
2652	pv[i].setKineticEnergyAndUpdate( ekin );
2653	}
2654	}
2655
2656	pvmx[4].setMass( incidentMass);
2657	pvmx[4].setMomentumAndUpdate( 0.0, 0.0, incidentTotalMomentum );
2658
2659	G4double ekin0 = pvmx[4].getKineticEnergy();
2660
2661	pvmx[5].setMass ( protonMass * targ);
2662	pvmx[5].setEnergy( protonMass * targ);
2663	pvmx[5].setKineticEnergy(0.);
2664	pvmx[5].setMomentum( 0.0, 0.0, 0.0 );
2665
2666	ekin = pvmx[4].getEnergy() + pvmx[5].getEnergy();
2667
2668	pvmx[6].Add( pvmx[4], pvmx[5] );
2669	pvmx[4].Lor( pvmx[4], pvmx[6] );
2670	pvmx[5].Lor( pvmx[5], pvmx[6] );
2671
2672	G4double tecm = pvmx[4].getEnergy() + pvmx[5].getEnergy();
2673
2674	pvmx[8].setZero();
2675
2676	G4double ekin1 = 0.0;
2677
2678	for( i=0; i < vecLen; i++ )
2679	{
2680	pvmx[8].Add( pvmx[8], pv[i] );
2681	ekin1 += pv[i].getKineticEnergy();
2682	ekin -= pv[i].getMass();
2683	}
2684
2685	if( vecLen > 1 && vecLen < 19 )
2686	{
2687	constantCrossSection = true;
2688	G4HEVector pw[18];
2689	for(i=0;i<vecLen;i++) pw[i] = pv[i];
2690	wgt = NBodyPhaseSpace( tecm, constantCrossSection, pw, vecLen );
2691	ekin = 0.0;
2692	for( i=0; i < vecLen; i++ )
2693	{
2694	pvmx[7].setMass( pw[i].getMass());
2695	pvmx[7].setMomentum( pw[i].getMomentum() );
2696	pvmx[7].SmulAndUpdate( pvmx[7], 1.);
2697	pvmx[7].Lor( pvmx[7], pvmx[5] );
2698	ekin += pvmx[7].getKineticEnergy();
2699	}
2700	teta = pvmx[8].Ang( pvmx[4] );
2701	if (verboseLevel > 1)
2702	G4cout << " vecLen > 1 && vecLen < 19 " << teta << " "
2703	<< ekin0 << " " << ekin1 << " " << ekin << G4endl;
2704	}
2705
2706	if( ekin1 != 0.0 )
2707	{
2708	pvmx[7].setZero();
2709	wgt = ekin/ekin1;
2710	ekin1 = 0.;
2711	for( i=0; i < vecLen; i++ )
2712	{
2713	pvMass = pv[i].getMass();
2714	ekin = pv[i].getKineticEnergy() * wgt;
2715	pv[i].setKineticEnergyAndUpdate( ekin );
2716	ekin1 += ekin;
2717	pvmx[7].Add( pvmx[7], pv[i] );
2718	}
2719	teta = pvmx[7].Ang( pvmx[4] );
2720	if (verboseLevel > 1)
2721	G4cout << " ekin1 != 0 " << teta << " " << ekin0 << " "
2722	<< ekin1 << G4endl;
2723	}
2724
2725	if(verboseLevel > 1)
2726	{
2727	G4cout << " After energy- correction " << G4endl;
2728	for(i=0; i<vecLen; i++) pv[i].Print(i);
2729	}
2730
2731	// Do some smearing in the transverse direction due to Fermi motion
2732
2733	G4double ry = G4UniformRand();
2734	G4double rz = G4UniformRand();
2735	G4double rx = twopi*rz;
2736	G4double a1 = std::sqrt(-2.0*std::log(ry));
2737	G4double rantarg1 = a1std::cos(rx)0.02*targ/G4double(vecLen);
2738	G4double rantarg2 = a1std::sin(rx)0.02*targ/G4double(vecLen);
2739
2740	for (i = 0; i < vecLen; i++)
2741	pv[i].setMomentum( pv[i].getMomentum().x()+rantarg1,
2742	pv[i].getMomentum().y()+rantarg2 );
2743
2744	if (verboseLevel > 1) {
2745	pvmx[7].setZero();
2746	for (i = 0; i < vecLen; i++) pvmx[7].Add( pvmx[7], pv[i] );
2747	teta = pvmx[7].Ang( pvmx[4] );
2748	G4cout << " After smearing " << teta << G4endl;
2749	}
2750
2751	// Rotate in the direction of the primary particle momentum (z-axis).
2752	// This does disturb our inclusive distributions somewhat, but it is
2753	// necessary for momentum conservation
2754
2755	// Also subtract binding energies and make some further corrections
2756	// if required
2757
2758	G4double dekin = 0.0;
2759	G4int npions = 0;
2760	G4double ek1 = 0.0;
2761	G4double alekw, xxh;
2762	G4double cfa = 0.025((atomicWeight-1.)/120.)std::exp(-(atomicWeight-1.)/120.);
2763	G4double alem[] = {1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00, 10.0};
2764	G4double val0[] = {0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65, 0.70};
2765
2766
2767	for (i = 0; i < vecLen; i++)
2768	{
2769	pv[i].Defs1( pv[i], pvI );
2770	if (atomicWeight > 1.5)
2771	{
2772	ekin = Amax( 1.e-6,pv[i].getKineticEnergy() - cfa( 1. + 0.5normal()));
2773	alekw = std::log( incidentKineticEnergy );
2774	xxh = 1.;
2775	xxh = 1.;
2776	if(incidentCode == pionPlusCode \|\| incidentCode == pionMinusCode)
2777	{
2778	if(pv[i].getCode() == pionZeroCode)
2779	{
2780	if(G4UniformRand() < std::log(atomicWeight))
2781	{
2782	if (alekw > alem[0])
2783	{
2784	for (j = 1; j < 8; j++)
2785	{
2786	if(alekw < alem[j]) break;
2787	}
2788	xxh = (val0[j]-val0[j-1])/(alem[j]-alem[j-1])*alekw
2789	+ val0[j-1] - (val0[j]-val0[j-1])/(alem[j]-alem[j-1])*alem[j-1];
2790	xxh = 1. - xxh;
2791	}
2792	}
2793	}
2794	}
2795	dekin += ekin*(1.-xxh);
2796	ekin *= xxh;
2797	pv[i].setKineticEnergyAndUpdate( ekin );
2798	pvCode = pv[i].getCode();
2799	if ((pvCode == pionPlusCode) \|\| (pvCode == pionMinusCode) \|\| (pvCode == pionZeroCode))
2800	{
2801	npions += 1;
2802	ek1 += ekin;
2803	}
2804	}
2805	}
2806	if( (ek1 > 0.0) && (npions > 0) )
2807	{
2808	dekin = 1.+dekin/ek1;
2809	for (i = 0; i < vecLen; i++)
2810	{
2811	pvCode = pv[i].getCode();
2812	if((pvCode == pionPlusCode) \|\| (pvCode == pionMinusCode) \|\| (pvCode == pionZeroCode))
2813	{
2814	ekin = Amax( 1.0e-6, pv[i].getKineticEnergy() * dekin );
2815	pv[i].setKineticEnergyAndUpdate( ekin );
2816	}
2817	}
2818	}
2819	if (verboseLevel > 1)
2820	{ G4cout << " Lab-System " << ek1 << " " << npions << G4endl;
2821	for (i=0; i<vecLen; i++) pv[i].Print(i);
2822	}
2823
2824	// Add black track particles
2825	// The total number of particles produced is restricted to 198
2826	// - this may have influence on very high energies
2827
2828	if (verboseLevel > 1)
2829	G4cout << " Evaporation " << atomicWeight << " " << excitationEnergyGNP
2830	<< " " << excitationEnergyDTA << G4endl;
2831
2832	G4double sprob = 0.;
2833	if (incidentKineticEnergy > 5. )
2834	// sprob = Amin( 1., (0.394-0.063std::log(atomicWeight))std::log(incidentKineticEnergy-4.) );
2835	sprob = Amin(1., 0.000314atomicWeightstd::log(incidentKineticEnergy-4.));
2836	if( atomicWeight > 1.5 && G4UniformRand() > sprob)
2837	{
2838
2839	G4double cost, sint, ekin2, ran, pp, eka;
2840	G4int spall(0), nbl(0);
2841
2842	// first add protons and neutrons
2843
2844	if( excitationEnergyGNP >= 0.001 )
2845	{
2846	// nbl = number of proton/neutron black track particles
2847	// tex is their total kinetic energy (GeV)
2848
2849	nbl = Poisson( (1.5+1.25targ)excitationEnergyGNP/
2850	(excitationEnergyGNP+excitationEnergyDTA));
2851	if( targ+nbl > atomicWeight ) nbl = (int)(atomicWeight - targ);
2852	if (verboseLevel > 1)
2853	G4cout << " evaporation " << targ << " " << nbl << " "
2854	<< sprob << G4endl;
2855	spall = targ;
2856	if( nbl > 0)
2857	{
2858	ekin = excitationEnergyGNP/nbl;
2859	ekin2 = 0.0;
2860	for( i=0; i<nbl; i++ )
2861	{
2862	if( G4UniformRand() < sprob ) continue;
2863	if( ekin2 > excitationEnergyGNP) break;
2864	ran = G4UniformRand();
2865	ekin1 = -ekinstd::log(ran) - cfa(1.0+0.5*normal());
2866	if (ekin1 < 0) ekin1 = -0.010*std::log(ran);
2867	ekin2 += ekin1;
2868	if( ekin2 > excitationEnergyGNP)
2869	ekin1 = Amax( 1.0e-6, excitationEnergyGNP-(ekin2-ekin1) );
2870	if( G4UniformRand() > (1.0-atomicNumber/(atomicWeight)))
2871	pv[vecLen].setDefinition( "Proton");
2872	else
2873	pv[vecLen].setDefinition( "Neutron" );
2874	spall++;
2875	cost = G4UniformRand() * 2.0 - 1.0;
2876	sint = std::sqrt(std::fabs(1.0-cost*cost));
2877	phi = twopi * G4UniformRand();
2878	pv[vecLen].setFlag( true ); // true is the same as IPA(i)<0
2879	pv[vecLen].setSide( -4 );
2880	pvMass = pv[vecLen].getMass();
2881	pv[vecLen].setTOF( 1.0 );
2882	pvEnergy = ekin1 + pvMass;
2883	pp = std::sqrt( std::fabs( pvEnergypvEnergy - pvMasspvMass ) );
2884	pv[vecLen].setMomentumAndUpdate( ppsintstd::sin(phi),
2885	ppsintstd::cos(phi),
2886	pp*cost );
2887	if (verboseLevel > 1) pv[vecLen].Print(vecLen);
2888	vecLen++;
2889	}
2890	if( (atomicWeight >= 10.0 ) && (incidentKineticEnergy <= 2.0) )
2891	{
2892	G4int ika, kk = 0;
2893	eka = incidentKineticEnergy;
2894	if( eka > 1.0 )eka *= eka;
2895	eka = Amax( 0.1, eka );
2896	ika = G4int(3.6std::exp((atomicNumberatomicNumber
2897	/atomicWeight-35.56)/6.45)/eka);
2898	if( ika > 0 )
2899	{
2900	for( i=(vecLen-1); i>=0; i-- )
2901	{
2902	if( (pv[i].getCode() == protonCode) && pv[i].getFlag() )
2903	{
2904	G4HEVector pTemp = pv[i];
2905	pv[i].setDefinition( "Neutron" );
2906	pv[i].setMomentumAndUpdate(pTemp.getMomentum());
2907	if (verboseLevel > 1) pv[i].Print(i);
2908	if( ++kk > ika ) break;
2909	}
2910	}
2911	}
2912	}
2913	}
2914	}
2915
2916	// Finished adding proton/neutron black track particles
2917	// now, try to add deuterons, tritons and alphas
2918
2919	if( excitationEnergyDTA >= 0.001 )
2920	{
2921	nbl = Poisson( (1.5+1.25targ)excitationEnergyDTA
2922	/(excitationEnergyGNP+excitationEnergyDTA));
2923
2924	// nbl is the number of deutrons, tritons, and alphas produced
2925
2926	if( nbl > 0 )
2927	{
2928	ekin = excitationEnergyDTA/nbl;
2929	ekin2 = 0.0;
2930	for( i=0; i<nbl; i++ )
2931	{
2932	if( G4UniformRand() < sprob ) continue;
2933	if( ekin2 > excitationEnergyDTA) break;
2934	ran = G4UniformRand();
2935	ekin1 = -ekinstd::log(ran)-cfa(1.+0.5*normal());
2936	if( ekin1 < 0.0 ) ekin1 = -0.010*std::log(ran);
2937	ekin2 += ekin1;
2938	if( ekin2 > excitationEnergyDTA )
2939	ekin1 = Amax( 1.0e-6, excitationEnergyDTA-(ekin2-ekin1));
2940	cost = G4UniformRand()*2.0 - 1.0;
2941	sint = std::sqrt(std::fabs(1.0-cost*cost));
2942	phi = twopi*G4UniformRand();
2943	ran = G4UniformRand();
2944	if( ran <= 0.60 )
2945	pv[vecLen].setDefinition( "Deuteron");
2946	else if (ran <= 0.90)
2947	pv[vecLen].setDefinition( "Triton" );
2948	else
2949	pv[vecLen].setDefinition( "Alpha" );
2950	spall += (int)(pv[vecLen].getMass() * 1.066);
2951	if( spall > atomicWeight ) break;
2952	pv[vecLen].setFlag( true ); // true is the same as IPA(i)<0
2953	pv[vecLen].setSide( -4 );
2954	pvMass = pv[vecLen].getMass();
2955	pv[vecLen].setTOF( 1.0 );
2956	pvEnergy = pvMass + ekin1;
2957	pp = std::sqrt( std::fabs( pvEnergypvEnergy - pvMasspvMass ) );
2958	pv[vecLen].setMomentumAndUpdate( ppsintstd::sin(phi),
2959	ppsintstd::cos(phi),
2960	pp*cost );
2961	if (verboseLevel > 1) pv[vecLen].Print(vecLen);
2962	vecLen++;
2963	}
2964	}
2965	}
2966	}
2967	if( centerOfMassEnergy <= (4.0+G4UniformRand()) )
2968	{
2969	for( i=0; i<vecLen; i++ )
2970	{
2971	G4double etb = pv[i].getKineticEnergy();
2972	if( etb >= incidentKineticEnergy )
2973	pv[i].setKineticEnergyAndUpdate( incidentKineticEnergy );
2974	}
2975	}
2976
2977	TuningOfHighEnergyCascading( pv, vecLen,
2978	incidentParticle, targetParticle,
2979	atomicWeight, atomicNumber);
2980
2981	// Calculate time delay for nuclear reactions
2982
2983	G4double tof = incidentTOF;
2984	if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0)
2985	&& (incidentKineticEnergy <= 0.2) )
2986	tof -= 500.0 * std::exp(-incidentKineticEnergy /0.04) * std::log( G4UniformRand() );
2987	for ( i=0; i < vecLen; i++)
2988	{
2989
2990	pv[i].setTOF ( tof );
2991	// vec[i].SetTOF ( tof );
2992	}
2993
2994	for(i=0; i<vecLen; i++)
2995	{
2996	if(pv[i].getName() == "KaonZero" \|\| pv[i].getName() == "AntiKaonZero")
2997	{
2998	pvmx[0] = pv[i];
2999	if(G4UniformRand() < 0.5) pv[i].setDefinition("KaonZeroShort");
3000	else pv[i].setDefinition("KaonZeroLong");
3001	pv[i].setMomentumAndUpdate(pvmx[0].getMomentum());
3002	}
3003	}
3004
3005	successful = true;
3006	delete [] pvmx;
3007	delete [] tempV;
3008	return;
3009	}
3010
3011	void
3012	G4HEInelastic::MediumEnergyCascading(G4bool &successful,
3013	G4HEVector pv[],
3014	G4int &vecLen,
3015	G4double &excitationEnergyGNP,
3016	G4double &excitationEnergyDTA,
3017	G4HEVector incidentParticle,
3018	G4HEVector targetParticle,
3019	G4double atomicWeight,
3020	G4double atomicNumber)
3021	{
3022	//
3023	// The multiplicity of particles produced in the first interaction has been
3024	// calculated in one of the FirstIntInNuc.... routines. The nuclear
3025	// cascading particles are parametrized from experimental data.
3026	// A simple single variable description E D3S/DP3= F(Q) with
3027	// Q^2 = (M*X)^2 + PT^2 is used. Final state kinematic is produced
3028	// by an FF-type iterative cascade method.
3029	// Nuclear evaporation particles are added at the end of the routine.
3030
3031	// All quantities on the G4HEVector Array pv are in GeV- units.
3032
3033	G4int protonCode = Proton.getCode();
3034	G4double protonMass = Proton.getMass();
3035	G4int neutronCode = Neutron.getCode();
3036	G4double kaonPlusMass = KaonPlus.getMass();
3037	G4int kaonPlusCode = KaonPlus.getCode();
3038	G4int kaonMinusCode = KaonMinus.getCode();
3039	G4int kaonZeroSCode = KaonZeroShort.getCode();
3040	G4int kaonZeroLCode = KaonZeroLong.getCode();
3041	G4int kaonZeroCode = KaonZero.getCode();
3042	G4int antiKaonZeroCode = AntiKaonZero.getCode();
3043	G4int pionPlusCode = PionPlus.getCode();
3044	G4int pionZeroCode = PionZero.getCode();
3045	G4int pionMinusCode = PionMinus.getCode();
3046	G4String mesonType = PionPlus.getType();
3047	G4String baryonType = Proton.getType();
3048	G4String antiBaryonType= AntiProton.getType();
3049
3050	G4double targetMass = targetParticle.getMass();
3051
3052	G4int incidentCode = incidentParticle.getCode();
3053	G4double incidentMass = incidentParticle.getMass();
3054	G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
3055	G4double incidentEnergy = incidentParticle.getEnergy();
3056	G4double incidentKineticEnergy = incidentParticle.getKineticEnergy();
3057	G4String incidentType = incidentParticle.getType();
3058	// G4double incidentTOF = incidentParticle.getTOF();
3059	G4double incidentTOF = 0.;
3060
3061	// some local variables
3062
3063	G4int i, j, l;
3064
3065	if(verboseLevel > 1)
3066	G4cout << " G4HEInelastic::MediumEnergyCascading " << G4endl;
3067
3068	// define annihilation channels.
3069
3070	G4bool annihilation = false;
3071	if (incidentCode < 0 && incidentType == antiBaryonType &&
3072	pv[0].getType() != antiBaryonType &&
3073	pv[1].getType() != antiBaryonType )
3074	{
3075	annihilation = true;
3076	}
3077
3078	successful = false;
3079
3080	G4double twsup[] = { 1., 1., 0.7, 0.5, 0.3, 0.2, 0.1, 0.0 };
3081
3082	if(annihilation) goto start;
3083	if(vecLen >= 8) goto start;
3084	if(incidentKineticEnergy < 1.) return;
3085	if( ( incidentCode == kaonPlusCode \|\| incidentCode == kaonMinusCode
3086	\|\| incidentCode == kaonZeroCode \|\| incidentCode == antiKaonZeroCode
3087	\|\| incidentCode == kaonZeroSCode \|\| incidentCode == kaonZeroLCode )
3088	&& ( G4UniformRand() < 0.5)) goto start;
3089	if(G4UniformRand() > twsup[vecLen-1]) goto start;
3090	return;
3091
3092	start:
3093
3094	if (annihilation)
3095	{ // do some corrections of incident particle kinematic
3096	G4double ekcor = Amax( 1., 1./incidentKineticEnergy);
3097	incidentKineticEnergy = 2targetMass + incidentKineticEnergy(1.+ekcor/atomicWeight);
3098	G4double excitation = NuclearExcitation(incidentKineticEnergy,
3099	atomicWeight,
3100	atomicNumber,
3101	excitationEnergyGNP,
3102	excitationEnergyDTA);
3103	incidentKineticEnergy -= excitation;
3104	if (incidentKineticEnergy < excitationEnergyDTA) incidentKineticEnergy = 0.;
3105	incidentEnergy = incidentKineticEnergy + incidentMass;
3106	incidentTotalMomentum =
3107	std::sqrt( Amax(0., incidentEnergyincidentEnergy - incidentMassincidentMass));
3108	}
3109
3110	G4HEVector pTemp;
3111	for(i = 2; i<vecLen; i++)
3112	{
3113	j = Imin(vecLen-1, (G4int)(2. + G4UniformRand()*(vecLen-2)));
3114	pTemp = pv[j];
3115	pv[j] = pv[i];
3116	pv[i] = pTemp;
3117	}
3118
3119	// randomize the first two leading particles
3120	// for kaon induced reactions only
3121	// (need from experimental data)
3122
3123	if( (incidentCode==kaonPlusCode \|\| incidentCode==kaonMinusCode \|\|
3124	incidentCode==kaonZeroCode \|\| incidentCode==antiKaonZeroCode \|\|
3125	incidentCode==kaonZeroSCode \|\| incidentCode==kaonZeroLCode)
3126	&& (G4UniformRand()>0.7) )
3127	{
3128	pTemp = pv[1];
3129	pv[1] = pv[0];
3130	pv[0] = pTemp;
3131	}
3132
3133	// mark leading particles for incident strange particles
3134	// and antibaryons, for all other we assume that the first
3135	// and second particle are the leading particles.
3136	// We need this later for kinematic aspects of strangeness
3137	// conservation.
3138
3139	G4int lead = 0;
3140	G4HEVector leadParticle;
3141	if( (incidentMass >= kaonPlusMass-0.05) && (incidentCode != protonCode)
3142	&& (incidentCode != neutronCode) )
3143	{
3144	G4double pMass = pv[0].getMass();
3145	G4int pCode = pv[0].getCode();
3146	if( (pMass >= kaonPlusMass-0.05) && (pCode != protonCode)
3147	&& (pCode != neutronCode) )
3148	{
3149	lead = pCode;
3150	leadParticle = pv[0];
3151	}
3152	else
3153	{
3154	pMass = pv[1].getMass();
3155	pCode = pv[1].getCode();
3156	if( (pMass >= kaonPlusMass-0.05) && (pCode != protonCode)
3157	&& (pCode != neutronCode) )
3158	{
3159	lead = pCode;
3160	leadParticle = pv[1];
3161	}
3162	}
3163	}
3164
3165	// Distribute particles in forward and backward hemispheres in center of
3166	// mass system. Incident particle goes in forward hemisphere.
3167
3168	G4HEVector pvI = incidentParticle; // for the incident particle
3169	pvI.setSide( 1 );
3170
3171	G4HEVector pvT = targetParticle; // for the target particle
3172	pvT.setMomentumAndUpdate( 0.0, 0.0, 0.0 );
3173	pvT.setSide( -1 );
3174	pvT.setTOF( -1.);
3175
3176
3177	G4double centerOfMassEnergy = std::sqrt( sqr(incidentMass)+sqr(targetMass)
3178	+2.0targetMassincidentEnergy );
3179	// G4double availableEnergy = centerOfMassEnergy - ( targetMass + incidentMass );
3180
3181	G4double tavai1 = centerOfMassEnergy/2.0 - incidentMass;
3182	G4double tavai2 = centerOfMassEnergy/2.0 - targetMass;
3183
3184	G4int ntb = 1;
3185	for( i=0; i < vecLen; i++ )
3186	{
3187	if (i == 0) pv[i].setSide( 1 );
3188	else if (i == 1) pv[i].setSide( -1 );
3189	else
3190	{ if( G4UniformRand() < 0.5 )
3191	{
3192	pv[i].setSide( -1 );
3193	ntb++;
3194	}
3195	else
3196	pv[i].setSide( 1 );
3197	}
3198	pv[i].setTOF( incidentTOF);
3199	}
3200	G4double tb = 2. * ntb;
3201	if (centerOfMassEnergy < (2. + G4UniformRand()))
3202	tb = (2. * ntb + vecLen)/2.;
3203
3204	if (verboseLevel > 1)
3205	{ G4cout << " pv Vector after Randomization " << vecLen << G4endl;
3206	pvI.Print(-1);
3207	pvT.Print(-1);
3208	for (i=0; i < vecLen ; i++) pv[i].Print(i);
3209	}
3210
3211	// Add particles from intranuclear cascade
3212	// nuclearCascadeCount = number of new secondaries
3213	// produced by nuclear cascading.
3214	// extraCount = number of nucleons within these new secondaries
3215
3216	G4double s, xtarg, ran;
3217	s = centerOfMassEnergy*centerOfMassEnergy;
3218	xtarg = Amax( 0.01, Amin( 0.75, 0.312 + 0.200 * std::log(std::log(s))
3219	+ std::pow(s,1.5)/6000.0 )
3220	(std::pow(atomicWeight, 0.33) - 1.0) tb);
3221
3222	G4int ntarg = Poisson( xtarg );
3223	G4int targ = 0;
3224
3225	if( ntarg > 0 )
3226	{
3227	G4double nucsup[] = { 1.00, 0.7, 0.5, 0.4, 0.35, 0.3 };
3228	G4double psup[] = { 3., 6., 20., 50., 100., 1000. };
3229	G4int momentumBin = 0;
3230	while( (momentumBin < 6) && (incidentTotalMomentum > psup[momentumBin]) )
3231	momentumBin++;
3232	momentumBin = Imin( 5, momentumBin );
3233
3234	// NOTE: in GENXPT, these new particles were given negative codes
3235	// here I use flag = true instead
3236
3237	for( i=0; i<ntarg; i++ )
3238	{
3239	if( G4UniformRand() < nucsup[momentumBin] )
3240	{
3241	if( G4UniformRand() > 1.0-atomicNumber/atomicWeight )
3242	pv[vecLen].setDefinition( "Proton" );
3243	else
3244	pv[vecLen].setDefinition( "Neutron" );
3245	targ++;
3246	}
3247	else
3248	{
3249	ran = G4UniformRand();
3250	if( ran < 0.33333 )
3251	pv[vecLen].setDefinition( "PionPlus");
3252	else if( ran < 0.66667 )
3253	pv[vecLen].setDefinition( "PionZero");
3254	else
3255	pv[vecLen].setDefinition( "PionMinus" );
3256	}
3257	pv[vecLen].setSide( -2 ); // backward cascade particles
3258	pv[vecLen].setFlag( true ); // true is the same as IPA(i)<0
3259	pv[vecLen].setTOF( incidentTOF );
3260	vecLen++;
3261	}
3262	}
3263
3264	// assume conservation of kinetic energy
3265	// in forward & backward hemispheres
3266
3267	G4int is, iskip;
3268	tavai1 = centerOfMassEnergy/2.;
3269	G4int iavai1 = 0;
3270
3271	for (i = 0; i < vecLen; i++)
3272	{
3273	if (pv[i].getSide() > 0)
3274	{
3275	tavai1 -= pv[i].getMass();
3276	iavai1++;
3277	}
3278	}
3279	if ( iavai1 == 0) return;
3280
3281	while( tavai1 <= 0.0 )
3282	{ // must eliminate a particle from the forward side
3283	iskip = G4int(G4UniformRand()*iavai1) + 1;
3284	is = 0;
3285	for( i=vecLen-1; i>=0; i-- )
3286	{
3287	if( pv[i].getSide() > 0 )
3288	{
3289	if (++is == iskip)
3290	{
3291	tavai1 += pv[i].getMass();
3292	iavai1--;
3293	if ( i != vecLen-1)
3294	{
3295	for( j=i; j<vecLen; j++ )
3296	{
3297	pv[j] = pv[j+1];
3298	}
3299	}
3300	if( --vecLen == 0 ) return; // all the secondaries except of the
3301	break; // --+
3302	} // \|
3303	} // v
3304	} // break goes down to here
3305	} // to the end of the for- loop.
3306
3307
3308	tavai2 = (targ+1)*centerOfMassEnergy/2.;
3309	G4int iavai2 = 0;
3310
3311	for (i = 0; i < vecLen; i++)
3312	{
3313	if (pv[i].getSide() < 0)
3314	{
3315	tavai2 -= pv[i].getMass();
3316	iavai2++;
3317	}
3318	}
3319	if (iavai2 == 0) return;
3320
3321	while( tavai2 <= 0.0 )
3322	{ // must eliminate a particle from the backward side
3323	iskip = G4int(G4UniformRand()*iavai2) + 1;
3324	is = 0;
3325	for( i = vecLen-1; i >= 0; i-- )
3326	{
3327	if( pv[i].getSide() < 0 )
3328	{
3329	if( ++is == iskip )
3330	{
3331	tavai2 += pv[i].getMass();
3332	iavai2--;
3333	if (pv[i].getSide() == -2) ntarg--;
3334	if (i != vecLen-1)
3335	{
3336	for( j=i; j<vecLen; j++)
3337	{
3338	pv[j] = pv[j+1];
3339	}
3340	}
3341	if (--vecLen == 0) return;
3342	break;
3343	}
3344	}
3345	}
3346	}
3347
3348	if (verboseLevel > 1) {
3349	G4cout << " pv Vector after Energy checks " << vecLen << " "
3350	<< tavai1 << " " << iavai1 << " " << tavai2 << " "
3351	<< iavai2 << " " << ntarg << G4endl;
3352	pvI.Print(-1);
3353	pvT.Print(-1);
3354	for (i=0; i < vecLen ; i++) pv[i].Print(i);
3355	}
3356
3357	// Define some vectors for Lorentz transformations
3358
3359	G4HEVector* pvmx = new G4HEVector [10];
3360
3361	pvmx[0].setMass( incidentMass );
3362	pvmx[0].setMomentumAndUpdate( 0.0, 0.0, incidentTotalMomentum );
3363	pvmx[1].setMass( protonMass);
3364	pvmx[1].setMomentumAndUpdate( 0.0, 0.0, 0.0 );
3365	pvmx[3].setMass( protonMass*(1+targ));
3366	pvmx[3].setMomentumAndUpdate( 0.0, 0.0, 0.0 );
3367	pvmx[4].setZero();
3368	pvmx[5].setZero();
3369	pvmx[7].setZero();
3370	pvmx[8].setZero();
3371	pvmx[8].setMomentum( 1.0, 0.0 );
3372	pvmx[2].Add( pvmx[0], pvmx[1] );
3373	pvmx[3].Add( pvmx[3], pvmx[0] );
3374	pvmx[0].Lor( pvmx[0], pvmx[2] );
3375	pvmx[1].Lor( pvmx[1], pvmx[2] );
3376
3377	if (verboseLevel > 1) {
3378	G4cout << " General Vectors after Definition " << G4endl;
3379	for (i=0; i<10; i++) pvmx[i].Print(i);
3380	}
3381
3382	// Main loop for 4-momentum generation - see Pitha-report (Aachen)
3383	// for a detailed description of the method.
3384	// Process the secondary particles in reverse order
3385
3386	G4double dndl[20];
3387	G4double binl[20];
3388	G4double pvMass, pvEnergy;
3389	G4int pvCode;
3390	G4double aspar, pt, phi, et, xval;
3391	G4double ekin = 0.;
3392	G4double ekin1 = 0.;
3393	G4double ekin2 = 0.;
3394	phi = G4UniformRand()*twopi;
3395	G4int npg = 0;
3396	G4int targ1 = 0; // No fragmentation model for nucleons
3397	for( i=vecLen-1; i>=0; i-- ) // from the intranuclear cascade. Mark
3398	{ // them with -3 and leave the loop.
3399	if( (pv[i].getSide() == -2) \|\| (i == 1) )
3400	{
3401	if ( pv[i].getType() == baryonType \|\|
3402	pv[i].getType() == antiBaryonType)
3403	{
3404	if( ++npg < 19 )
3405	{
3406	pv[i].setSide( -3 );
3407	targ1++;
3408	continue; // leave the for loop !!
3409	}
3410	}
3411	}
3412
3413	// Set pt and phi values - they are changed somewhat in the
3414	// iteration loop.
3415	// Set mass parameter for lambda fragmentation model
3416
3417	G4double maspar[] = { 0.75, 0.70, 0.65, 0.60, 0.50, 0.40, 0.75, 0.20};
3418	G4double bp[] = { 3.50, 3.50, 3.50, 6.00, 5.00, 4.00, 3.50, 3.50};
3419	G4double ptex[] = { 1.70, 1.70, 1.50, 1.70, 1.40, 1.20, 1.70, 1.20};
3420	// Set parameters for lambda simulation:
3421	// pt is the average transverse momentum
3422	// aspar the is average transverse mass
3423
3424	pvMass = pv[i].getMass();
3425	j = 2;
3426	if ( pv[i].getType() == mesonType ) j = 1;
3427	if ( pv[i].getMass() < 0.4 ) j = 0;
3428	if ( i <= 1 ) j += 3;
3429	if (pv[i].getSide() <= -2) j = 6;
3430	if (j == 6 && (pv[i].getType() == baryonType \|\| pv[i].getType()==antiBaryonType) ) j = 7;
3431	pt = Amax(0.001, std::sqrt(std::pow(-std::log(1.-G4UniformRand())/bp[j],ptex[j])));
3432	aspar = maspar[j];
3433	phi = G4UniformRand()*twopi;
3434	pv[i].setMomentum( ptstd::cos(phi), ptstd::sin(phi) ); // set x- and y-momentum
3435
3436	for( j=0; j<20; j++ ) binl[j] = j/(19.*pt); // set the lambda - bins.
3437
3438	if( pv[i].getSide() > 0 )
3439	et = pvmx[0].getEnergy();
3440	else
3441	et = pvmx[1].getEnergy();
3442
3443	dndl[0] = 0.0;
3444
3445	// Start of outer iteration loop
3446
3447	G4int outerCounter = 0, innerCounter = 0; // three times.
3448	G4bool eliminateThisParticle = true;
3449	G4bool resetEnergies = true;
3450	while( ++outerCounter < 3 )
3451	{
3452	for( l=1; l<20; l++ )
3453	{
3454	xval = (binl[l]+binl[l-1])/2.; // x = lambda /GeV
3455	if( xval > 1./pt )
3456	dndl[l] = dndl[l-1];
3457	else
3458	dndl[l] = dndl[l-1] +
3459	aspar/std::sqrt( std::pow((1.+asparxvalasparxval),3) )
3460	(binl[l]-binl[l-1]) * et /
3461	std::sqrt( ptxvaletptxvalet + ptpt + pvMass*pvMass );
3462	}
3463
3464	// Start of inner iteration loop
3465
3466	innerCounter = 0; // try this not more than 7 times.
3467	while( ++innerCounter < 7 )
3468	{
3469	l = 1;
3470	ran = G4UniformRand()*dndl[19];
3471	while( ( ran >= dndl[l] ) && ( l < 20 ) )l++;
3472	l = Imin( 19, l );
3473	xval = Amin( 1.0, pt(binl[l-1] + G4UniformRand()(binl[l]-binl[l-1]) ) );
3474	if( pv[i].getSide() < 0 ) xval *= -1.;
3475	pv[i].setMomentumAndUpdate( xval*et ); // set the z-momentum
3476	pvEnergy = pv[i].getEnergy();
3477	if( pv[i].getSide() > 0 ) // forward side
3478	{
3479	if ( i < 2 )
3480	{
3481	ekin = tavai1 - ekin1;
3482	if (ekin < 0.) ekin = 0.04*std::fabs(normal());
3483	G4double pp1 = pv[i].Length();
3484	if (pp1 >= 1.e-6)
3485	{
3486	G4double pp = std::sqrt(ekin(ekin+2pvMass));
3487	pp = Amax(0.,pppp-ptpt);
3488	pp = std::sqrt(pp)*pv[i].getSide()/std::fabs(G4double(pv[i].getSide()));
3489	pv[i].setMomentumAndUpdate( pp );
3490	}
3491	else
3492	{
3493	pv[i].setMomentum(0.,0.,0.);
3494	pv[i].setKineticEnergyAndUpdate( ekin);
3495	}
3496	pvmx[4].Add( pvmx[4], pv[i]);
3497	outerCounter = 2;
3498	resetEnergies = false;
3499	eliminateThisParticle = false;
3500	break;
3501	}
3502	else if( (ekin1+pvEnergy-pvMass) < 0.95*tavai1 )
3503	{
3504	pvmx[4].Add( pvmx[4], pv[i] );
3505	ekin1 += pvEnergy - pvMass;
3506	pvmx[6].Add( pvmx[4], pvmx[5] );
3507	pvmx[6].setMomentum( 0.0 );
3508	outerCounter = 2; // leave outer loop
3509	eliminateThisParticle = false; // don't eliminate this particle
3510	resetEnergies = false;
3511	break; // next particle
3512	}
3513	if( innerCounter > 5 ) break; // leave inner loop
3514
3515	if( tavai2 >= pvMass )
3516	{ // switch sides
3517	pv[i].setSide( -1 );
3518	tavai1 += pvMass;
3519	tavai2 -= pvMass;
3520	iavai2++;
3521	}
3522	}
3523	else
3524	{ // backward side
3525	xval = Amin(0.999,0.95+0.05*targ/20.0);
3526	if( (ekin2+pvEnergy-pvMass) < xval*tavai2 )
3527	{
3528	pvmx[5].Add( pvmx[5], pv[i] );
3529	ekin2 += pvEnergy - pvMass;
3530	pvmx[6].Add( pvmx[4], pvmx[5] );
3531	pvmx[6].setMomentum( 0.0 ); // set z-momentum
3532	outerCounter = 2; // leave outer iteration
3533	eliminateThisParticle = false; // don't eliminate this particle
3534	resetEnergies = false;
3535	break; // leave inner iteration
3536	}
3537	if( innerCounter > 5 )break; // leave inner iteration
3538
3539	if( tavai1 >= pvMass )
3540	{ // switch sides
3541	pv[i].setSide( 1 );
3542	tavai1 -= pvMass;
3543	tavai2 += pvMass;
3544	iavai2--;
3545	}
3546	}
3547	pv[i].setMomentum( pv[i].getMomentum().x() * 0.9,
3548	pv[i].getMomentum().y() * 0.9);
3549	pt *= 0.9;
3550	dndl[19] *= 0.9;
3551	} // closes inner loop
3552
3553	if (resetEnergies)
3554	{
3555	ekin1 = 0.0;
3556	ekin2 = 0.0;
3557	pvmx[4].setZero();
3558	pvmx[5].setZero();
3559	if (verboseLevel > 1)
3560	G4cout << " Reset energies for index " << i << G4endl;
3561	for( l=i+1; l < vecLen; l++ )
3562	{
3563	if( (pv[l].getMass() < protonMass) \|\| (pv[l].getSide() > 0) )
3564	{
3565	pvEnergy = Amax( pv[l].getMass(), 0.95*pv[l].getEnergy()
3566	+ 0.05*pv[l].getMass() );
3567	pv[l].setEnergyAndUpdate( pvEnergy );
3568	if( pv[l].getSide() > 0)
3569	{
3570	ekin1 += pv[l].getKineticEnergy();
3571	pvmx[4].Add( pvmx[4], pv[l] );
3572	}
3573	else
3574	{
3575	ekin2 += pv[l].getKineticEnergy();
3576	pvmx[5].Add( pvmx[5], pv[l] );
3577	}
3578	}
3579	}
3580	}
3581	} // closes outer iteration
3582
3583	if( eliminateThisParticle ) // not enough energy,
3584	{ // eliminate this particle
3585	if (verboseLevel > 1)
3586	{
3587	G4cout << " Eliminate particle with index " << i << G4endl;
3588	pv[i].Print(i);
3589	}
3590	for( j=i; j < vecLen; j++ )
3591	{ // shift down
3592	pv[j] = pv[j+1];
3593	}
3594	vecLen--;
3595	if(vecLen < 2) return;
3596	i++;
3597	pvmx[6].Add( pvmx[4], pvmx[5] );
3598	pvmx[6].setMomentum( 0.0 ); // set z-momentum
3599	}
3600	} // closes main for loop
3601	if (verboseLevel > 1)
3602	{ G4cout << " pv Vector after lambda fragmentation " << vecLen << G4endl;
3603	pvI.Print(-1);
3604	pvT.Print(-1);
3605	for (i=0; i < vecLen ; i++) pv[i].Print(i);
3606	for (i=0; i < 10; i++) pvmx[i].Print(i);
3607	}
3608
3609	// Backward nucleons produced with a cluster model
3610
3611	pvmx[6].Lor( pvmx[3], pvmx[2] );
3612	pvmx[6].Sub( pvmx[6], pvmx[4] );
3613	pvmx[6].Sub( pvmx[6], pvmx[5] );
3614	if (verboseLevel > 1) pvmx[6].Print(6);
3615
3616	npg = 0;
3617	G4double rmb0 = 0.;
3618	G4double rmb;
3619	G4double wgt;
3620	G4bool constantCrossSection = true;
3621	for (i = 0; i < vecLen; i++)
3622	{
3623	if(pv[i].getSide() == -3)
3624	{
3625	npg++;
3626	rmb0 += pv[i].getMass();
3627	}
3628	}
3629	if( targ1 == 1 \|\| npg < 2)
3630	{ // target particle is the only backward nucleon
3631	ekin = Amin( tavai2-ekin2, centerOfMassEnergy/2.0-protonMass );
3632	if( ekin < 0.04 ) ekin = 0.04 * std::fabs( normal() );
3633	G4double pp = std::sqrt(ekin(ekin+2pv[1].getMass()));
3634	G4double pp1 = pvmx[6].Length();
3635	if(pp1 < 1.e-6)
3636	{
3637	pv[1].setKineticEnergyAndUpdate(ekin);
3638	}
3639	else
3640	{
3641	pv[1].setMomentum(pvmx[6].getMomentum());
3642	pv[1].SmulAndUpdate(pv[1],pp/pp1);
3643	}
3644	pvmx[5].Add( pvmx[5], pv[1] );
3645	}
3646	else
3647	{
3648	G4double cpar[] = { 0.6, 0.6, 0.35, 0.15, 0.10 };
3649	G4double gpar[] = { 2.6, 2.6, 1.80, 1.30, 1.20 };
3650
3651	G4int tempCount = Imin( 5, targ1 ) - 1;
3652
3653	rmb = rmb0 + std::pow(-std::log(1.0-G4UniformRand()), cpar[tempCount])/gpar[tempCount];
3654	pvEnergy = pvmx[6].getEnergy();
3655	if ( rmb > pvEnergy ) rmb = pvEnergy;
3656	pvmx[6].setMass( rmb );
3657	pvmx[6].setEnergyAndUpdate( pvEnergy );
3658	pvmx[6].Smul( pvmx[6], -1. );
3659	if (verboseLevel > 1) {
3660	G4cout << " General Vectors before input to NBodyPhaseSpace "
3661	<< targ1 << " " << tempCount << " " << rmb0 << " "
3662	<< rmb << " " << pvEnergy << G4endl;
3663	for (i=0; i<10; i++) pvmx[i].Print(i);
3664	}
3665
3666	// tempV contains the backward nucleons
3667
3668	G4HEVector* tempV = new G4HEVector[18];
3669	npg = 0;
3670	for( i=0; i < vecLen; i++ )
3671	{
3672	if( pv[i].getSide() == -3 ) tempV[npg++] = pv[i];
3673	}
3674
3675	wgt = NBodyPhaseSpace( pvmx[6].getMass(), constantCrossSection, tempV, npg );
3676
3677	npg = 0;
3678	for( i=0; i < vecLen; i++ )
3679	{
3680	if( pv[i].getSide() == -3 )
3681	{
3682	pv[i].setMomentum( tempV[npg++].getMomentum());
3683	pv[i].SmulAndUpdate( pv[i], 1.);
3684	pv[i].Lor( pv[i], pvmx[6] );
3685	pvmx[5].Add( pvmx[5], pv[i] );
3686	}
3687	}
3688	delete [] tempV;
3689	}
3690	if( vecLen <= 2 )
3691	{
3692	successful = false;
3693	return;
3694	}
3695
3696	// Lorentz transformation in lab system
3697
3698	targ = 0;
3699	for( i=0; i < vecLen; i++ )
3700	{
3701	if( pv[i].getType() == baryonType )targ++;
3702	if( pv[i].getType() == antiBaryonType )targ++;
3703	pv[i].Lor( pv[i], pvmx[1] );
3704	}
3705	targ = Imax( 1, targ );
3706
3707	G4bool dum(0);
3708	if( lead )
3709	{
3710	for( i=0; i<vecLen; i++ )
3711	{
3712	if( pv[i].getCode() == lead )
3713	{
3714	dum = false;
3715	break;
3716	}
3717	}
3718	if( dum )
3719	{
3720	i = 0;
3721
3722	if( ( (leadParticle.getType() == baryonType \|\|
3723	leadParticle.getType() == antiBaryonType)
3724	&& (pv[1].getType() == baryonType \|\|
3725	pv[1].getType() == antiBaryonType))
3726	\|\| ( (leadParticle.getType() == mesonType)
3727	&& (pv[1].getType() == mesonType)))
3728	{
3729	i = 1;
3730	}
3731	ekin = pv[i].getKineticEnergy();
3732	pv[i] = leadParticle;
3733	if( pv[i].getFlag() )
3734	pv[i].setTOF( -1.0 );
3735	else
3736	pv[i].setTOF( 1.0 );
3737	pv[i].setKineticEnergyAndUpdate( ekin );
3738	}
3739	}
3740
3741	pvmx[3].setMass( incidentMass);
3742	pvmx[3].setMomentumAndUpdate( 0.0, 0.0, incidentTotalMomentum );
3743
3744	G4double ekin0 = pvmx[3].getKineticEnergy();
3745
3746	pvmx[4].setMass ( protonMass * targ);
3747	pvmx[4].setEnergy( protonMass * targ);
3748	pvmx[4].setMomentum(0.,0.,0.);
3749	pvmx[4].setKineticEnergy(0.);
3750
3751	ekin = pvmx[3].getEnergy() + pvmx[4].getEnergy();
3752
3753	pvmx[5].Add( pvmx[3], pvmx[4] );
3754	pvmx[3].Lor( pvmx[3], pvmx[5] );
3755	pvmx[4].Lor( pvmx[4], pvmx[5] );
3756
3757	G4double tecm = pvmx[3].getEnergy() + pvmx[4].getEnergy();
3758
3759	pvmx[7].setZero();
3760
3761	ekin1 = 0.0;
3762	G4double teta;
3763
3764	for( i=0; i < vecLen; i++ )
3765	{
3766	pvmx[7].Add( pvmx[7], pv[i] );
3767	ekin1 += pv[i].getKineticEnergy();
3768	ekin -= pv[i].getMass();
3769	}
3770
3771	if( vecLen > 1 && vecLen < 19 )
3772	{
3773	constantCrossSection = true;
3774	G4HEVector pw[18];
3775	for(i=0;i<vecLen;i++) pw[i] = pv[i];
3776	wgt = NBodyPhaseSpace( tecm, constantCrossSection, pw, vecLen );
3777	ekin = 0.0;
3778	for( i=0; i < vecLen; i++ )
3779	{
3780	pvmx[6].setMass( pw[i].getMass());
3781	pvmx[6].setMomentum( pw[i].getMomentum() );
3782	pvmx[6].SmulAndUpdate( pvmx[6], 1.);
3783	pvmx[6].Lor( pvmx[6], pvmx[4] );
3784	ekin += pvmx[6].getKineticEnergy();
3785	}
3786	teta = pvmx[7].Ang( pvmx[3] );
3787	if (verboseLevel > 1)
3788	G4cout << " vecLen > 1 && vecLen < 19 " << teta << " " << ekin0
3789	<< " " << ekin1 << " " << ekin << G4endl;
3790	}
3791
3792	if( ekin1 != 0.0 )
3793	{
3794	pvmx[6].setZero();
3795	wgt = ekin/ekin1;
3796	ekin1 = 0.;
3797	for( i=0; i < vecLen; i++ )
3798	{
3799	pvMass = pv[i].getMass();
3800	ekin = pv[i].getKineticEnergy() * wgt;
3801	pv[i].setKineticEnergyAndUpdate( ekin );
3802	ekin1 += ekin;
3803	pvmx[6].Add( pvmx[6], pv[i] );
3804	}
3805	teta = pvmx[6].Ang( pvmx[3] );
3806	if (verboseLevel > 1)
3807	G4cout << " ekin1 != 0 " << teta << " " << ekin0 << " "
3808	<< ekin1 << G4endl;
3809	}
3810
3811	// Do some smearing in the transverse direction due to Fermi motion.
3812
3813	G4double ry = G4UniformRand();
3814	G4double rz = G4UniformRand();
3815	G4double rx = twopi*rz;
3816	G4double a1 = std::sqrt(-2.0*std::log(ry));
3817	G4double rantarg1 = a1std::cos(rx)0.02*targ/G4double(vecLen);
3818	G4double rantarg2 = a1std::sin(rx)0.02*targ/G4double(vecLen);
3819
3820	for (i = 0; i < vecLen; i++)
3821	pv[i].setMomentum( pv[i].getMomentum().x()+rantarg1,
3822	pv[i].getMomentum().y()+rantarg2 );
3823
3824	if (verboseLevel > 1) {
3825	pvmx[6].setZero();
3826	for (i = 0; i < vecLen; i++) pvmx[6].Add( pvmx[6], pv[i] );
3827	teta = pvmx[6].Ang( pvmx[3] );
3828	G4cout << " After smearing " << teta << G4endl;
3829	}
3830
3831	// Rotate in the direction of the primary particle momentum (z-axis).
3832	// This does disturb our inclusive distributions somewhat, but it is
3833	// necessary for momentum conservation.
3834
3835	// Also subtract binding energies and make some further corrections
3836	// if required.
3837
3838	G4double dekin = 0.0;
3839	G4int npions = 0;
3840	G4double ek1 = 0.0;
3841	G4double alekw, xxh;
3842	G4double cfa = 0.025((atomicWeight-1.)/120.)std::exp(-(atomicWeight-1.)/120.);
3843	G4double alem[] = {1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00, 10.00};
3844	G4double val0[] = {0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65, 0.70};
3845
3846
3847	for (i = 0; i < vecLen; i++)
3848	{
3849	pv[i].Defs1( pv[i], pvI );
3850	if (atomicWeight > 1.5)
3851	{
3852	ekin = Amax( 1.e-6,pv[i].getKineticEnergy() - cfa( 1. + 0.5normal()));
3853	alekw = std::log( incidentKineticEnergy );
3854	xxh = 1.;
3855	if(incidentCode == pionPlusCode \|\| incidentCode == pionMinusCode)
3856	{
3857	if(pv[i].getCode() == pionZeroCode)
3858	{
3859	if(G4UniformRand() < std::log(atomicWeight))
3860	{
3861	if (alekw > alem[0])
3862	{
3863	for (j = 1; j < 8; j++)
3864	{
3865	if(alekw < alem[j]) break;
3866	}
3867	xxh = (val0[j]-val0[j-1])/(alem[j]-alem[j-1])*alekw
3868	+ val0[j-1] - (val0[j]-val0[j-1])/(alem[j]-alem[j-1])*alem[j-1];
3869	xxh = 1. - xxh;
3870	}
3871	}
3872	}
3873	}
3874	dekin += ekin*(1.-xxh);
3875	ekin *= xxh;
3876	pv[i].setKineticEnergyAndUpdate( ekin );
3877	pvCode = pv[i].getCode();
3878	if ((pvCode == pionPlusCode) \|\| (pvCode == pionMinusCode) \|\| (pvCode == pionZeroCode))
3879	{
3880	npions += 1;
3881	ek1 += ekin;
3882	}
3883	}
3884	}
3885	if( (ek1 > 0.0) && (npions > 0) )
3886	{
3887	dekin = 1.+dekin/ek1;
3888	for (i = 0; i < vecLen; i++)
3889	{
3890	pvCode = pv[i].getCode();
3891	if((pvCode == pionPlusCode) \|\| (pvCode == pionMinusCode) \|\| (pvCode == pionZeroCode))
3892	{
3893	ekin = Amax( 1.0e-6, pv[i].getKineticEnergy() * dekin );
3894	pv[i].setKineticEnergyAndUpdate( ekin );
3895	}
3896	}
3897	}
3898	if (verboseLevel > 1)
3899	{ G4cout << " Lab-System " << ek1 << " " << npions << G4endl;
3900	for (i=0; i<vecLen; i++) pv[i].Print(i);
3901	}
3902
3903	// Add black track particles
3904	// The total number of particles produced is restricted to 198
3905	// this may have influence on very high energies
3906
3907	if (verboseLevel > 1) G4cout << " Evaporation " << atomicWeight << " " <<
3908	excitationEnergyGNP << " " << excitationEnergyDTA << G4endl;
3909
3910	if( atomicWeight > 1.5 )
3911	{
3912
3913	G4double sprob, cost, sint, pp, eka;
3914	G4int spall(0), nbl(0);
3915	// sprob is the probability of self-absorption in heavy molecules
3916
3917	if( incidentKineticEnergy < 5.0 )
3918	sprob = 0.0;
3919	else
3920	// sprob = Amin( 1.0, 0.6*std::log(incidentKineticEnergy-4.0) );
3921	sprob = Amin(1., 0.000314atomicWeightstd::log(incidentKineticEnergy-4.));
3922
3923	// First add protons and neutrons
3924
3925	if( excitationEnergyGNP >= 0.001 )
3926	{
3927	// nbl = number of proton/neutron black track particles
3928	// tex is their total kinetic energy (GeV)
3929
3930	nbl = Poisson( (1.5+1.25targ)excitationEnergyGNP/
3931	(excitationEnergyGNP+excitationEnergyDTA));
3932	if( targ+nbl > atomicWeight ) nbl = (int)(atomicWeight - targ);
3933	if (verboseLevel > 1)
3934	G4cout << " evaporation " << targ << " " << nbl << " "
3935	<< sprob << G4endl;
3936	spall = targ;
3937	if( nbl > 0)
3938	{
3939	ekin = excitationEnergyGNP/nbl;
3940	ekin2 = 0.0;
3941	for( i=0; i<nbl; i++ )
3942	{
3943	if( G4UniformRand() < sprob ) continue;
3944	if( ekin2 > excitationEnergyGNP) break;
3945	ran = G4UniformRand();
3946	ekin1 = -ekinstd::log(ran) - cfa(1.0+0.5*normal());
3947	if (ekin1 < 0) ekin1 = -0.010*std::log(ran);
3948	ekin2 += ekin1;
3949	if( ekin2 > excitationEnergyGNP)
3950	ekin1 = Amax( 1.0e-6, excitationEnergyGNP-(ekin2-ekin1) );
3951	if( G4UniformRand() > (1.0-atomicNumber/(atomicWeight)))
3952	pv[vecLen].setDefinition( "Proton");
3953	else
3954	pv[vecLen].setDefinition( "Neutron");
3955	spall++;
3956	cost = G4UniformRand() * 2.0 - 1.0;
3957	sint = std::sqrt(std::fabs(1.0-cost*cost));
3958	phi = twopi * G4UniformRand();
3959	pv[vecLen].setFlag( true ); // true is the same as IPA(i)<0
3960	pv[vecLen].setSide( -4 );
3961	pvMass = pv[vecLen].getMass();
3962	pv[vecLen].setTOF( 1.0 );
3963	pvEnergy = ekin1 + pvMass;
3964	pp = std::sqrt( std::fabs( pvEnergypvEnergy - pvMasspvMass ) );
3965	pv[vecLen].setMomentumAndUpdate( ppsintstd::sin(phi),
3966	ppsintstd::cos(phi),
3967	pp*cost );
3968	if (verboseLevel > 1) pv[vecLen].Print(vecLen);
3969	vecLen++;
3970	}
3971	if( (atomicWeight >= 10.0 ) && (incidentKineticEnergy <= 2.0) )
3972	{
3973	G4int ika, kk = 0;
3974	eka = incidentKineticEnergy;
3975	if( eka > 1.0 )eka *= eka;
3976	eka = Amax( 0.1, eka );
3977	ika = G4int(3.6std::exp((atomicNumberatomicNumber
3978	/atomicWeight-35.56)/6.45)/eka);
3979	if( ika > 0 )
3980	{
3981	for( i=(vecLen-1); i>=0; i-- )
3982	{
3983	if( (pv[i].getCode() == protonCode) && pv[i].getFlag() )
3984	{
3985	pTemp = pv[i];
3986	pv[i].setDefinition( "Neutron");
3987	pv[i].setMomentumAndUpdate(pTemp.getMomentum());
3988	if (verboseLevel > 1) pv[i].Print(i);
3989	if( ++kk > ika ) break;
3990	}
3991	}
3992	}
3993	}
3994	}
3995	}
3996
3997	// Finished adding proton/neutron black track particles
3998	// now, try to add deuterons, tritons and alphas
3999
4000	if( excitationEnergyDTA >= 0.001 )
4001	{
4002	nbl = Poisson( (1.5+1.25targ)excitationEnergyDTA
4003	/(excitationEnergyGNP+excitationEnergyDTA));
4004
4005	// nbl is the number of deutrons, tritons, and alphas produced
4006
4007	if( nbl > 0 )
4008	{
4009	ekin = excitationEnergyDTA/nbl;
4010	ekin2 = 0.0;
4011	for( i=0; i<nbl; i++ )
4012	{
4013	if( G4UniformRand() < sprob ) continue;
4014	if( ekin2 > excitationEnergyDTA) break;
4015	ran = G4UniformRand();
4016	ekin1 = -ekinstd::log(ran)-cfa(1.+0.5*normal());
4017	if( ekin1 < 0.0 ) ekin1 = -0.010*std::log(ran);
4018	ekin2 += ekin1;
4019	if( ekin2 > excitationEnergyDTA)
4020	ekin1 = Amax( 1.0e-6, excitationEnergyDTA-(ekin2-ekin1));
4021	cost = G4UniformRand()*2.0 - 1.0;
4022	sint = std::sqrt(std::fabs(1.0-cost*cost));
4023	phi = twopi*G4UniformRand();
4024	ran = G4UniformRand();
4025	if( ran <= 0.60 )
4026	pv[vecLen].setDefinition( "Deuteron");
4027	else if (ran <= 0.90)
4028	pv[vecLen].setDefinition( "Triton");
4029	else
4030	pv[vecLen].setDefinition( "Alpha");
4031	spall += (int)(pv[vecLen].getMass() * 1.066);
4032	if( spall > atomicWeight ) break;
4033	pv[vecLen].setFlag( true ); // true is the same as IPA(i)<0
4034	pv[vecLen].setSide( -4 );
4035	pvMass = pv[vecLen].getMass();
4036	pv[vecLen].setSide( pv[vecLen].getCode());
4037	pv[vecLen].setTOF( 1.0 );
4038	pvEnergy = pvMass + ekin1;
4039	pp = std::sqrt( std::fabs( pvEnergypvEnergy - pvMasspvMass ) );
4040	pv[vecLen].setMomentumAndUpdate( ppsintstd::sin(phi),
4041	ppsintstd::cos(phi),
4042	pp*cost );
4043	if (verboseLevel > 1) pv[vecLen].Print(vecLen);
4044	vecLen++;
4045	}
4046	}
4047	}
4048	}
4049	if( centerOfMassEnergy <= (4.0+G4UniformRand()) )
4050	{
4051	for( i=0; i<vecLen; i++ )
4052	{
4053	G4double etb = pv[i].getKineticEnergy();
4054	if( etb >= incidentKineticEnergy )
4055	pv[i].setKineticEnergyAndUpdate( incidentKineticEnergy );
4056	}
4057	}
4058
4059	// Calculate time delay for nuclear reactions
4060
4061	G4double tof = incidentTOF;
4062	if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0)
4063	&& (incidentKineticEnergy <= 0.2) )
4064	tof -= 500.0 * std::exp(-incidentKineticEnergy /0.04) * std::log( G4UniformRand() );
4065	for ( i=0; i < vecLen; i++)
4066	{
4067
4068	pv[i].setTOF ( tof );
4069	// vec[i].SetTOF ( tof );
4070	}
4071
4072	for(i=0; i<vecLen; i++)
4073	{
4074	if(pv[i].getName() == "KaonZero" \|\| pv[i].getName() == "AntiKaonZero")
4075	{
4076	pvmx[0] = pv[i];
4077	if(G4UniformRand() < 0.5) pv[i].setDefinition("KaonZeroShort");
4078	else pv[i].setDefinition("KaonZeroLong");
4079	pv[i].setMomentumAndUpdate(pvmx[0].getMomentum());
4080	}
4081	}
4082
4083	successful = true;
4084	delete [] pvmx;
4085	return;
4086	}
4087
4088	void
4089	G4HEInelastic::MediumEnergyClusterProduction(G4bool &successful,
4090	G4HEVector pv[],
4091	G4int &vecLen,
4092	G4double &excitationEnergyGNP,
4093	G4double &excitationEnergyDTA,
4094	G4HEVector incidentParticle,
4095	G4HEVector targetParticle,
4096	G4double atomicWeight,
4097	G4double atomicNumber)
4098	{
4099	// For low multiplicity in the first intranuclear interaction the cascading
4100	// process as described in G4HEInelastic::MediumEnergyCascading does not work
4101	// satisfactorily. From experimental data it is strongly suggested to use
4102	// a two- body resonance model.
4103	//
4104	// All quantities on the G4HEVector Array pv are in GeV- units.
4105
4106	G4int protonCode = Proton.getCode();
4107	G4double protonMass = Proton.getMass();
4108	G4int neutronCode = Neutron.getCode();
4109	G4double kaonPlusMass = KaonPlus.getMass();
4110	G4int pionPlusCode = PionPlus.getCode();
4111	G4int pionZeroCode = PionZero.getCode();
4112	G4int pionMinusCode = PionMinus.getCode();
4113	G4String mesonType = PionPlus.getType();
4114	G4String baryonType = Proton.getType();
4115	G4String antiBaryonType= AntiProton.getType();
4116
4117	G4double targetMass = targetParticle.getMass();
4118
4119	G4int incidentCode = incidentParticle.getCode();
4120	G4double incidentMass = incidentParticle.getMass();
4121	G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
4122	G4double incidentEnergy = incidentParticle.getEnergy();
4123	G4double incidentKineticEnergy = incidentParticle.getKineticEnergy();
4124	G4String incidentType = incidentParticle.getType();
4125	// G4double incidentTOF = incidentParticle.getTOF();
4126	G4double incidentTOF = 0.;
4127
4128	// some local variables
4129
4130	G4int i, j;
4131
4132	if(verboseLevel > 1) G4cout << " G4HEInelastic::MediumEnergyClusterProduction " << G4endl;
4133
4134	if (incidentTotalMomentum < 0.01)
4135	{
4136	successful = false;
4137	return;
4138	}
4139	G4double centerOfMassEnergy = std::sqrt( sqr(incidentMass) + sqr(targetMass)
4140	+2.targetMassincidentEnergy);
4141
4142	G4HEVector pvI = incidentParticle; // for the incident particle
4143	pvI.setSide( 1 );
4144
4145	G4HEVector pvT = targetParticle; // for the target particle
4146	pvT.setMomentumAndUpdate( 0.0, 0.0, 0.0 );
4147	pvT.setSide( -1 );
4148	pvT.setTOF( -1.);
4149
4150	// Distribute particles in forward and backward hemispheres. Note that
4151	// only low multiplicity events from FirstIntInNuc.... should go into
4152	// this routine.
4153
4154	G4int targ = 0;
4155	G4int ifor = 0;
4156	G4int iback = 0;
4157	G4int pvCode;
4158	G4double pvMass, pvEnergy;
4159
4160	pv[0].setSide( 1 );
4161	pv[1].setSide( -1 );
4162	for(i = 0; i < vecLen; i++)
4163	{
4164	if (i > 1)
4165	{
4166	if( G4UniformRand() < 0.5)
4167	{
4168	pv[i].setSide( 1 );
4169	if (++ifor > 18)
4170	{
4171	pv[i].setSide( -1 );
4172	ifor--;
4173	iback++;
4174	}
4175	}
4176	else
4177	{
4178	pv[i].setSide( -1 );
4179	if (++iback > 18)
4180	{
4181	pv[i].setSide( 1 );
4182	ifor++;
4183	iback--;
4184	}
4185	}
4186	}
4187
4188	pvCode = pv[i].getCode();
4189
4190	if ( ( (incidentCode == protonCode) \|\| (incidentCode == neutronCode)
4191	\|\| (incidentType == mesonType) )
4192	&& ( (pvCode == pionPlusCode) \|\| (pvCode == pionMinusCode) )
4193	&& ( (G4UniformRand() < (10.-incidentTotalMomentum)/6.) )
4194	&& ( (G4UniformRand() < atomicWeight/300.) ) )
4195	{
4196	if (G4UniformRand() > atomicNumber/atomicWeight)
4197	pv[i].setDefinition( "Neutron");
4198	else
4199	pv[i].setDefinition( "Proton");
4200	targ++;
4201	}
4202	pv[i].setTOF( incidentTOF );
4203	}
4204	G4double tb = 2. * iback;
4205	if (centerOfMassEnergy < (2+G4UniformRand())) tb = (2.*iback + vecLen)/2.;
4206
4207	G4double nucsup[] = { 1.0, 0.8, 0.6, 0.5, 0.4};
4208
4209	G4double xtarg = Amax(0.01, (0.312+0.2std::log(std::log(centerOfMassEnergycenterOfMassEnergy)))
4210	* (std::pow(atomicWeight,0.33)-1.) * tb);
4211	G4int ntarg = Poisson(xtarg);
4212	if (ntarg > 0)
4213	{
4214	G4int ipx = Imin(4, (G4int)(incidentTotalMomentum/3.));
4215	for (i=0; i < ntarg; i++)
4216	{
4217	if (G4UniformRand() < nucsup[ipx] )
4218	{
4219	if (G4UniformRand() < (1.- atomicNumber/atomicWeight))
4220	pv[vecLen].setDefinition( "Neutron");
4221	else
4222	pv[vecLen].setDefinition( "Proton");
4223	targ++;
4224	}
4225	else
4226	{
4227	G4double ran = G4UniformRand();
4228	if (ran < 0.3333 )
4229	pv[vecLen].setDefinition( "PionPlus");
4230	else if (ran < 0.6666)
4231	pv[vecLen].setDefinition( "PionZero");
4232	else
4233	pv[vecLen].setDefinition( "PionMinus");
4234	}
4235	pv[vecLen].setSide( -2 );
4236	pv[vecLen].setFlag( true );
4237	pv[vecLen].setTOF( incidentTOF );
4238	vecLen++;
4239	}
4240	}
4241
4242	// Mark leading particles for incident strange particles and antibaryons,
4243	// for all other we assume that the first and second particle are the
4244	// leading particles.
4245	// We need this later for kinematic aspects of strangeness conservation.
4246
4247	G4int lead = 0;
4248	G4HEVector leadParticle;
4249	if( (incidentMass >= kaonPlusMass-0.05) && (incidentCode != protonCode)
4250	&& (incidentCode != neutronCode) )
4251	{
4252	G4double pMass = pv[0].getMass();
4253	G4int pCode = pv[0].getCode();
4254	if( (pMass >= kaonPlusMass-0.05) && (pCode != protonCode)
4255	&& (pCode != neutronCode) )
4256	{
4257	lead = pCode;
4258	leadParticle = pv[0];
4259	}
4260	else
4261	{
4262	pMass = pv[1].getMass();
4263	pCode = pv[1].getCode();
4264	if( (pMass >= kaonPlusMass-0.05) && (pCode != protonCode)
4265	&& (pCode != neutronCode) )
4266	{
4267	lead = pCode;
4268	leadParticle = pv[1];
4269	}
4270	}
4271	}
4272
4273	if (verboseLevel > 1) {
4274	G4cout << " pv Vector after initialization " << vecLen << G4endl;
4275	pvI.Print(-1);
4276	pvT.Print(-1);
4277	for (i=0; i < vecLen ; i++) pv[i].Print(i);
4278	}
4279
4280	G4double tavai = 0.;
4281	for(i=0;i<vecLen;i++) if(pv[i].getSide() != -2) tavai += pv[i].getMass();
4282
4283	while (tavai > centerOfMassEnergy)
4284	{
4285	for (i=vecLen-1; i >= 0; i--)
4286	{
4287	if (pv[i].getSide() != -2)
4288	{
4289	tavai -= pv[i].getMass();
4290	if( i != vecLen-1)
4291	{
4292	for (j=i; j < vecLen; j++)
4293	{
4294	pv[j] = pv[j+1];
4295	}
4296	}
4297	if ( --vecLen < 2)
4298	{
4299	successful = false;
4300	return;
4301	}
4302	break;
4303	}
4304	}
4305	}
4306
4307	// Now produce 3 Clusters:
4308	// 1. forward cluster
4309	// 2. backward meson cluster
4310	// 3. backward nucleon cluster
4311
4312	G4double rmc0 = 0., rmd0 = 0., rme0 = 0.;
4313	G4int ntc = 0, ntd = 0, nte = 0;
4314
4315	for (i=0; i < vecLen; i++)
4316	{
4317	if(pv[i].getSide() > 0)
4318	{
4319	if(ntc < 17)
4320	{
4321	rmc0 += pv[i].getMass();
4322	ntc++;
4323	}
4324	else
4325	{
4326	if(ntd < 17)
4327	{
4328	pv[i].setSide(-1);
4329	rmd0 += pv[i].getMass();
4330	ntd++;
4331	}
4332	else
4333	{
4334	pv[i].setSide(-2);
4335	rme0 += pv[i].getMass();
4336	nte++;
4337	}
4338	}
4339	}
4340	else if (pv[i].getSide() == -1)
4341	{
4342	if(ntd < 17)
4343	{
4344	rmd0 += pv[i].getMass();
4345	ntd++;
4346	}
4347	else
4348	{
4349	pv[i].setSide(-2);
4350	rme0 += pv[i].getMass();
4351	nte++;
4352	}
4353	}
4354	else
4355	{
4356	rme0 += pv[i].getMass();
4357	nte++;
4358	}
4359	}
4360
4361	G4double cpar[] = {0.6, 0.6, 0.35, 0.15, 0.10};
4362	G4double gpar[] = {2.6, 2.6, 1.80, 1.30, 1.20};
4363
4364	G4double rmc = rmc0, rmd = rmd0, rme = rme0;
4365	G4int ntc1 = Imin(4,ntc-1);
4366	G4int ntd1 = Imin(4,ntd-1);
4367	G4int nte1 = Imin(4,nte-1);
4368	if (ntc > 1) rmc = rmc0 + std::pow(-std::log(1.-G4UniformRand()),cpar[ntc1])/gpar[ntc1];
4369	if (ntd > 1) rmd = rmd0 + std::pow(-std::log(1.-G4UniformRand()),cpar[ntd1])/gpar[ntd1];
4370	if (nte > 1) rme = rme0 + std::pow(-std::log(1.-G4UniformRand()),cpar[nte1])/gpar[nte1];
4371	while( (rmc+rmd) > centerOfMassEnergy)
4372	{
4373	if ((rmc == rmc0) && (rmd == rmd0))
4374	{
4375	rmd = 0.999centerOfMassEnergy/(rmc+rmd);
4376	rmc = 0.999centerOfMassEnergy/(rmc+rmd);
4377	}
4378	else
4379	{
4380	rmc = 0.1rmc0 + 0.9rmc;
4381	rmd = 0.1rmd0 + 0.9rmd;
4382	}
4383	}
4384	if(verboseLevel > 1)
4385	G4cout << " Cluster Masses: " << ntc << " " << rmc << " " << ntd << " "
4386	<< rmd << " " << nte << " " << rme << G4endl;
4387
4388
4389	G4HEVector* pvmx = new G4HEVector[11];
4390
4391	pvmx[1].setMass( incidentMass);
4392	pvmx[1].setMomentumAndUpdate( 0., 0., incidentTotalMomentum);
4393	pvmx[2].setMass( targetMass);
4394	pvmx[2].setMomentumAndUpdate( 0., 0., 0.);
4395	pvmx[0].Add( pvmx[1], pvmx[2] );
4396	pvmx[1].Lor( pvmx[1], pvmx[0] );
4397	pvmx[2].Lor( pvmx[2], pvmx[0] );
4398
4399	G4double pf = std::sqrt(Amax(0.0001, sqr(sqr(centerOfMassEnergy) + rmdrmd -rmcrmc)
4400	- 4sqr(centerOfMassEnergy)rmdrmd))/(2.centerOfMassEnergy);
4401	pvmx[3].setMass( rmc );
4402	pvmx[4].setMass( rmd );
4403	pvmx[3].setEnergy( std::sqrt(pfpf + rmcrmc) );
4404	pvmx[4].setEnergy( std::sqrt(pfpf + rmdrmd) );
4405
4406	G4double tvalue = -MAXFLOAT;
4407	G4double bvalue = Amax(0.01, 4.0 + 1.6*std::log(incidentTotalMomentum));
4408	if (bvalue != 0.0) tvalue = std::log(G4UniformRand())/bvalue;
4409	G4double pin = pvmx[1].Length();
4410	G4double tacmin = sqr( pvmx[1].getEnergy() - pvmx[3].getEnergy()) - sqr( pin - pf);
4411	G4double ctet = Amax(-1., Amin(1., 1.+2.(tvalue-tacmin)/Amax(1.e-10, 4.pin*pf)));
4412	G4double stet = std::sqrt(Amax(0., 1.0 - ctet*ctet));
4413	G4double phi = twopi * G4UniformRand();
4414	pvmx[3].setMomentum( pf * stet * std::sin(phi),
4415	pf * stet * std::cos(phi),
4416	pf * ctet );
4417	pvmx[4].Smul( pvmx[3], -1.);
4418
4419	if (nte > 0)
4420	{
4421	G4double ekit1 = 0.04;
4422	G4double ekit2 = 0.6;
4423	G4double gaval = 1.2;
4424	if (incidentKineticEnergy <= 5.)
4425	{
4426	ekit1 *= sqr(incidentKineticEnergy)/25.;
4427	ekit2 *= sqr(incidentKineticEnergy)/25.;
4428	}
4429	G4double avalue = (1.-gaval)/(std::pow(ekit2, 1.-gaval)-std::pow(ekit1, 1.-gaval));
4430	for (i=0; i < vecLen; i++)
4431	{
4432	if (pv[i].getSide() == -2)
4433	{
4434	G4double ekit = std::pow(G4UniformRand()*(1.-gaval)/avalue +std::pow(ekit1, 1.-gaval),
4435	1./(1.-gaval));
4436	pv[i].setKineticEnergyAndUpdate( ekit );
4437	ctet = Amax(-1., Amin(1., std::log(2.23*G4UniformRand()+0.383)/0.96));
4438	stet = std::sqrt( Amax( 0.0, 1. - ctet*ctet ));
4439	phi = G4UniformRand()*twopi;
4440	G4double pp = pv[i].Length();
4441	pv[i].setMomentum( pp * stet * std::sin(phi),
4442	pp * stet * std::cos(phi),
4443	pp * ctet );
4444	pv[i].Lor( pv[i], pvmx[0] );
4445	}
4446	}
4447	}
4448	// pvmx[1] = pvmx[3];
4449	// pvmx[2] = pvmx[4];
4450	pvmx[5].SmulAndUpdate( pvmx[3], -1.);
4451	pvmx[6].SmulAndUpdate( pvmx[4], -1.);
4452
4453	if (verboseLevel > 1) {
4454	G4cout << " General vectors before Phase space Generation " << G4endl;
4455	for (i=0; i<7; i++) pvmx[i].Print(i);
4456	}
4457
4458
4459	G4HEVector* tempV = new G4HEVector[18];
4460	G4bool constantCrossSection = true;
4461	G4double wgt;
4462	G4int npg;
4463
4464	if (ntc > 1)
4465	{
4466	npg = 0;
4467	for (i=0; i < vecLen; i++)
4468	{
4469	if (pv[i].getSide() > 0)
4470	{
4471	tempV[npg++] = pv[i];
4472	if(verboseLevel > 1) pv[i].Print(i);
4473	}
4474	}
4475	wgt = NBodyPhaseSpace( pvmx[3].getMass(), constantCrossSection, tempV, npg);
4476
4477	npg = 0;
4478	for (i=0; i < vecLen; i++)
4479	{
4480	if (pv[i].getSide() > 0)
4481	{
4482	pv[i].setMomentum( tempV[npg++].getMomentum());
4483	pv[i].SmulAndUpdate( pv[i], 1. );
4484	pv[i].Lor( pv[i], pvmx[5] );
4485	if(verboseLevel > 1) pv[i].Print(i);
4486	}
4487	}
4488	}
4489	else if(ntc == 1)
4490	{
4491	for(i=0; i<vecLen; i++)
4492	{
4493	if(pv[i].getSide() > 0) pv[i].setMomentumAndUpdate(pvmx[3].getMomentum());
4494	if(verboseLevel > 1) pv[i].Print(i);
4495	}
4496	}
4497	else
4498	{
4499	}
4500
4501	if (ntd > 1)
4502	{
4503	npg = 0;
4504	for (i=0; i < vecLen; i++)
4505	{
4506	if (pv[i].getSide() == -1)
4507	{
4508	tempV[npg++] = pv[i];
4509	if(verboseLevel > 1) pv[i].Print(i);
4510	}
4511	}
4512	wgt = NBodyPhaseSpace( pvmx[4].getMass(), constantCrossSection, tempV, npg);
4513
4514	npg = 0;
4515	for (i=0; i < vecLen; i++)
4516	{
4517	if (pv[i].getSide() == -1)
4518	{
4519	pv[i].setMomentum( tempV[npg++].getMomentum());
4520	pv[i].SmulAndUpdate( pv[i], 1.);
4521	pv[i].Lor( pv[i], pvmx[6] );
4522	if(verboseLevel > 1) pv[i].Print(i);
4523	}
4524	}
4525	}
4526	else if(ntd == 1)
4527	{
4528	for(i=0; i<vecLen; i++)
4529	{
4530	if(pv[i].getSide() == -1) pv[i].setMomentumAndUpdate(pvmx[4].getMomentum());
4531	if(verboseLevel > 1) pv[i].Print(i);
4532	}
4533	}
4534	else
4535	{
4536	}
4537
4538	if(verboseLevel > 1)
4539	{
4540	G4cout << " Vectors after PhaseSpace generation " << G4endl;
4541	for(i=0;i<vecLen; i++) pv[i].Print(i);
4542	}
4543
4544	// Lorentz transformation in lab system
4545
4546	targ = 0;
4547	for( i=0; i < vecLen; i++ )
4548	{
4549	if( pv[i].getType() == baryonType )targ++;
4550	if( pv[i].getType() == antiBaryonType )targ++;
4551	pv[i].Lor( pv[i], pvmx[2] );
4552	}
4553	if (targ <1) targ =1;
4554
4555	if(verboseLevel > 1) {
4556	G4cout << " Transformation in Lab- System " << G4endl;
4557	for(i=0; i<vecLen; i++) pv[i].Print(i);
4558	}
4559
4560	G4bool dum(0);
4561	G4double ekin, teta;
4562
4563	if( lead )
4564	{
4565	for( i=0; i<vecLen; i++ )
4566	{
4567	if( pv[i].getCode() == lead )
4568	{
4569	dum = false;
4570	break;
4571	}
4572	}
4573	if( dum )
4574	{
4575	i = 0;
4576
4577	if( ( (leadParticle.getType() == baryonType \|\|
4578	leadParticle.getType() == antiBaryonType)
4579	&& (pv[1].getType() == baryonType \|\|
4580	pv[1].getType() == antiBaryonType))
4581	\|\| ( (leadParticle.getType() == mesonType)
4582	&& (pv[1].getType() == mesonType)))
4583	{
4584	i = 1;
4585	}
4586
4587	ekin = pv[i].getKineticEnergy();
4588	pv[i] = leadParticle;
4589	if( pv[i].getFlag() )
4590	pv[i].setTOF( -1.0 );
4591	else
4592	pv[i].setTOF( 1.0 );
4593	pv[i].setKineticEnergyAndUpdate( ekin );
4594	}
4595	}
4596
4597	pvmx[4].setMass( incidentMass);
4598	pvmx[4].setMomentumAndUpdate( 0.0, 0.0, incidentTotalMomentum );
4599
4600	G4double ekin0 = pvmx[4].getKineticEnergy();
4601
4602	pvmx[5].setMass ( protonMass * targ);
4603	pvmx[5].setMomentumAndUpdate( 0.0, 0.0, 0.0 );
4604
4605	ekin = pvmx[4].getEnergy() + pvmx[5].getEnergy();
4606
4607	pvmx[6].Add( pvmx[4], pvmx[5] );
4608	pvmx[4].Lor( pvmx[4], pvmx[6] );
4609	pvmx[5].Lor( pvmx[5], pvmx[6] );
4610
4611	G4double tecm = pvmx[4].getEnergy() + pvmx[5].getEnergy();
4612
4613	pvmx[8].setZero();
4614
4615	G4double ekin1 = 0.0;
4616
4617	for( i=0; i < vecLen; i++ )
4618	{
4619	pvmx[8].Add( pvmx[8], pv[i] );
4620	ekin1 += pv[i].getKineticEnergy();
4621	ekin -= pv[i].getMass();
4622	}
4623
4624	if( vecLen > 1 && vecLen < 19 )
4625	{
4626	constantCrossSection = true;
4627	G4HEVector pw[18];
4628	for(i=0;i<vecLen;i++) pw[i] = pv[i];
4629	wgt = NBodyPhaseSpace( tecm, constantCrossSection, pw, vecLen );
4630	ekin = 0.0;
4631	for( i=0; i < vecLen; i++ )
4632	{
4633	pvmx[7].setMass( pw[i].getMass());
4634	pvmx[7].setMomentum( pw[i].getMomentum() );
4635	pvmx[7].SmulAndUpdate( pvmx[7], 1.);
4636	pvmx[7].Lor( pvmx[7], pvmx[5] );
4637	ekin += pvmx[7].getKineticEnergy();
4638	}
4639	teta = pvmx[8].Ang( pvmx[4] );
4640	if (verboseLevel > 1)
4641	G4cout << " vecLen > 1 && vecLen < 19 " << teta << " " << ekin0
4642	<< " " << ekin1 << " " << ekin << G4endl;
4643	}
4644
4645	if( ekin1 != 0.0 )
4646	{
4647	pvmx[7].setZero();
4648	wgt = ekin/ekin1;
4649	ekin1 = 0.;
4650	for( i=0; i < vecLen; i++ )
4651	{
4652	pvMass = pv[i].getMass();
4653	ekin = pv[i].getKineticEnergy() * wgt;
4654	pv[i].setKineticEnergyAndUpdate( ekin );
4655	ekin1 += ekin;
4656	pvmx[7].Add( pvmx[7], pv[i] );
4657	}
4658	teta = pvmx[7].Ang( pvmx[4] );
4659	if (verboseLevel > 1)
4660	G4cout << " ekin1 != 0 " << teta << " " << ekin0 << " "
4661	<< ekin1 << G4endl;
4662	}
4663
4664	// Do some smearing in the transverse direction due to Fermi motion.
4665
4666	G4double ry = G4UniformRand();
4667	G4double rz = G4UniformRand();
4668	G4double rx = twopi*rz;
4669	G4double a1 = std::sqrt(-2.0*std::log(ry));
4670	G4double rantarg1 = a1std::cos(rx)0.02*targ/G4double(vecLen);
4671	G4double rantarg2 = a1std::sin(rx)0.02*targ/G4double(vecLen);
4672
4673	for (i = 0; i < vecLen; i++)
4674	pv[i].setMomentum( pv[i].getMomentum().x()+rantarg1,
4675	pv[i].getMomentum().y()+rantarg2 );
4676
4677	if (verboseLevel > 1) {
4678	pvmx[7].setZero();
4679	for (i = 0; i < vecLen; i++) pvmx[7].Add( pvmx[7], pv[i] );
4680	teta = pvmx[7].Ang( pvmx[4] );
4681	G4cout << " After smearing " << teta << G4endl;
4682	}
4683
4684	// Rotate in the direction of the primary particle momentum (z-axis).
4685	// This does disturb our inclusive distributions somewhat, but it is
4686	// necessary for momentum conservation.
4687
4688	// Also subtract binding energies and make some further corrections
4689	// if required.
4690
4691	G4double dekin = 0.0;
4692	G4int npions = 0;
4693	G4double ek1 = 0.0;
4694	G4double alekw, xxh;
4695	G4double cfa = 0.025((atomicWeight-1.)/120.)std::exp(-(atomicWeight-1.)/120.);
4696	G4double alem[] = {1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00};
4697	G4double val0[] = {0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65};
4698
4699
4700	for (i = 0; i < vecLen; i++)
4701	{
4702	pv[i].Defs1( pv[i], pvI );
4703	if (atomicWeight > 1.5)
4704	{
4705	ekin = Amax( 1.e-6,pv[i].getKineticEnergy() - cfa( 1. + 0.5normal()));
4706	alekw = std::log( incidentKineticEnergy );
4707	xxh = 1.;
4708	xxh = 1.;
4709	if(incidentCode == pionPlusCode \|\| incidentCode == pionMinusCode)
4710	{
4711	if(pv[i].getCode() == pionZeroCode)
4712	{
4713	if(G4UniformRand() < std::log(atomicWeight))
4714	{
4715	if (alekw > alem[0])
4716	{
4717	for (j = 1; j < 8; j++)
4718	{
4719	if(alekw < alem[j]) break;
4720	}
4721	xxh = (val0[j]-val0[j-1])/(alem[j]-alem[j-1])*alekw
4722	+ val0[j-1] - (val0[j]-val0[j-1])/(alem[j]-alem[j-1])*alem[j-1];
4723	xxh = 1. - xxh;
4724	}
4725	}
4726	}
4727	}
4728	dekin += ekin*(1.-xxh);
4729	ekin *= xxh;
4730	pv[i].setKineticEnergyAndUpdate( ekin );
4731	pvCode = pv[i].getCode();
4732	if ((pvCode == pionPlusCode) \|\| (pvCode == pionMinusCode) \|\| (pvCode == pionZeroCode))
4733	{
4734	npions += 1;
4735	ek1 += ekin;
4736	}
4737	}
4738	}
4739	if( (ek1 > 0.0) && (npions > 0) )
4740	{
4741	dekin = 1.+dekin/ek1;
4742	for (i = 0; i < vecLen; i++)
4743	{
4744	pvCode = pv[i].getCode();
4745	if((pvCode == pionPlusCode) \|\| (pvCode == pionMinusCode) \|\| (pvCode == pionZeroCode))
4746	{
4747	ekin = Amax( 1.0e-6, pv[i].getKineticEnergy() * dekin );
4748	pv[i].setKineticEnergyAndUpdate( ekin );
4749	}
4750	}
4751	}
4752	if (verboseLevel > 1)
4753	{ G4cout << " Lab-System " << ek1 << " " << npions << G4endl;
4754	for (i=0; i<vecLen; i++) pv[i].Print(i);
4755	}
4756
4757	// Add black track particles
4758	// The total number of particles produced is restricted to 198
4759	// this may have influence on very high energies
4760
4761	if (verboseLevel > 1)
4762	G4cout << " Evaporation " << atomicWeight << " "
4763	<< excitationEnergyGNP << " " << excitationEnergyDTA << G4endl;
4764
4765	if( atomicWeight > 1.5 )
4766	{
4767
4768	G4double sprob, cost, sint, ekin2, ran, pp, eka;
4769	G4int spall(0), nbl(0);
4770	// sprob is the probability of self-absorption in heavy molecules
4771
4772	if( incidentKineticEnergy < 5.0 )
4773	sprob = 0.0;
4774	else
4775	// sprob = Amin( 1.0, 0.6*std::log(incidentKineticEnergy-4.0) );
4776	sprob = Amin(1., 0.000314atomicWeightstd::log(incidentKineticEnergy-4.));
4777	// First add protons and neutrons
4778
4779	if( excitationEnergyGNP >= 0.001 )
4780	{
4781	// nbl = number of proton/neutron black track particles
4782	// tex is their total kinetic energy (GeV)
4783
4784	nbl = Poisson( (1.5+1.25targ)excitationEnergyGNP/
4785	(excitationEnergyGNP+excitationEnergyDTA));
4786	if( targ+nbl > atomicWeight ) nbl = (int)(atomicWeight - targ);
4787	if (verboseLevel > 1)
4788	G4cout << " evaporation " << targ << " " << nbl << " "
4789	<< sprob << G4endl;
4790	spall = targ;
4791	if( nbl > 0)
4792	{
4793	ekin = excitationEnergyGNP/nbl;
4794	ekin2 = 0.0;
4795	for( i=0; i<nbl; i++ )
4796	{
4797	if( G4UniformRand() < sprob ) continue;
4798	if( ekin2 > excitationEnergyGNP) break;
4799	ran = G4UniformRand();
4800	ekin1 = -ekinstd::log(ran) - cfa(1.0+0.5*normal());
4801	if (ekin1 < 0) ekin1 = -0.010*std::log(ran);
4802	ekin2 += ekin1;
4803	if( ekin2 > excitationEnergyGNP )
4804	ekin1 = Amax( 1.0e-6, excitationEnergyGNP-(ekin2-ekin1) );
4805	if( G4UniformRand() > (1.0-atomicNumber/(atomicWeight)))
4806	pv[vecLen].setDefinition( "Proton");
4807	else
4808	pv[vecLen].setDefinition( "Neutron");
4809	spall++;
4810	cost = G4UniformRand() * 2.0 - 1.0;
4811	sint = std::sqrt(std::fabs(1.0-cost*cost));
4812	phi = twopi * G4UniformRand();
4813	pv[vecLen].setFlag( true ); // true is the same as IPA(i)<0
4814	pv[vecLen].setSide( -4 );
4815	pvMass = pv[vecLen].getMass();
4816	pv[vecLen].setTOF( 1.0 );
4817	pvEnergy = ekin1 + pvMass;
4818	pp = std::sqrt( std::fabs( pvEnergypvEnergy - pvMasspvMass ) );
4819	pv[vecLen].setMomentumAndUpdate( ppsintstd::sin(phi),
4820	ppsintstd::cos(phi),
4821	pp*cost );
4822	if (verboseLevel > 1) pv[vecLen].Print(vecLen);
4823	vecLen++;
4824	}
4825	if( (atomicWeight >= 10.0 ) && (incidentKineticEnergy <= 2.0) )
4826	{
4827	G4int ika, kk = 0;
4828	eka = incidentKineticEnergy;
4829	if( eka > 1.0 )eka *= eka;
4830	eka = Amax( 0.1, eka );
4831	ika = G4int(3.6std::exp((atomicNumberatomicNumber
4832	/atomicWeight-35.56)/6.45)/eka);
4833	if( ika > 0 )
4834	{
4835	for( i=(vecLen-1); i>=0; i-- )
4836	{
4837	if( (pv[i].getCode() == protonCode) && pv[i].getFlag() )
4838	{
4839	G4HEVector pTemp = pv[i];
4840	pv[i].setDefinition( "Neutron");
4841	pv[i].setMomentumAndUpdate(pTemp.getMomentum());
4842	if (verboseLevel > 1) pv[i].Print(i);
4843	if( ++kk > ika ) break;
4844	}
4845	}
4846	}
4847	}
4848	}
4849	}
4850
4851	// Finished adding proton/neutron black track particles
4852	// now, try to add deuterons, tritons and alphas
4853
4854	if( excitationEnergyDTA >= 0.001 )
4855	{
4856	nbl = Poisson( (1.5+1.25targ)excitationEnergyDTA
4857	/(excitationEnergyGNP+excitationEnergyDTA));
4858
4859	// nbl is the number of deutrons, tritons, and alphas produced
4860
4861	if( nbl > 0 )
4862	{
4863	ekin = excitationEnergyDTA/nbl;
4864	ekin2 = 0.0;
4865	for( i=0; i<nbl; i++ )
4866	{
4867	if( G4UniformRand() < sprob ) continue;
4868	if( ekin2 > excitationEnergyDTA) break;
4869	ran = G4UniformRand();
4870	ekin1 = -ekinstd::log(ran)-cfa(1.+0.5*normal());
4871	if( ekin1 < 0.0 ) ekin1 = -0.010*std::log(ran);
4872	ekin2 += ekin1;
4873	if( ekin2 > excitationEnergyDTA)
4874	ekin1 = Amax( 1.0e-6, excitationEnergyDTA-(ekin2-ekin1));
4875	cost = G4UniformRand()*2.0 - 1.0;
4876	sint = std::sqrt(std::fabs(1.0-cost*cost));
4877	phi = twopi*G4UniformRand();
4878	ran = G4UniformRand();
4879	if( ran <= 0.60 )
4880	pv[vecLen].setDefinition( "Deuteron");
4881	else if (ran <= 0.90)
4882	pv[vecLen].setDefinition( "Triton");
4883	else
4884	pv[vecLen].setDefinition( "Alpha");
4885	spall += (int)(pv[vecLen].getMass() * 1.066);
4886	if( spall > atomicWeight ) break;
4887	pv[vecLen].setFlag( true ); // true is the same as IPA(i)<0
4888	pv[vecLen].setSide( -4 );
4889	pvMass = pv[vecLen].getMass();
4890	pv[vecLen].setTOF( 1.0 );
4891	pvEnergy = pvMass + ekin1;
4892	pp = std::sqrt( std::fabs( pvEnergypvEnergy - pvMasspvMass ) );
4893	pv[vecLen].setMomentumAndUpdate( ppsintstd::sin(phi),
4894	ppsintstd::cos(phi),
4895	pp*cost );
4896	if (verboseLevel > 1) pv[vecLen].Print(vecLen);
4897	vecLen++;
4898	}
4899	}
4900	}
4901	}
4902	if( centerOfMassEnergy <= (4.0+G4UniformRand()) )
4903	{
4904	for( i=0; i<vecLen; i++ )
4905	{
4906	G4double etb = pv[i].getKineticEnergy();
4907	if( etb >= incidentKineticEnergy )
4908	pv[i].setKineticEnergyAndUpdate( incidentKineticEnergy );
4909	}
4910	}
4911
4912	// Calculate time delay for nuclear reactions
4913
4914	G4double tof = incidentTOF;
4915	if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0)
4916	&& (incidentKineticEnergy <= 0.2) )
4917	tof -= 500.0 * std::exp(-incidentKineticEnergy /0.04) * std::log( G4UniformRand() );
4918	for ( i=0; i < vecLen; i++)
4919	{
4920
4921	pv[i].setTOF ( tof );
4922	// vec[i].SetTOF ( tof );
4923	}
4924
4925	for(i=0; i<vecLen; i++)
4926	{
4927	if(pv[i].getName() == "KaonZero" \|\| pv[i].getName() == "AntiKaonZero")
4928	{
4929	pvmx[0] = pv[i];
4930	if(G4UniformRand() < 0.5) pv[i].setDefinition("KaonZeroShort");
4931	else pv[i].setDefinition("KaonZeroLong");
4932	pv[i].setMomentumAndUpdate(pvmx[0].getMomentum());
4933	}
4934	}
4935
4936	successful = true;
4937	delete [] pvmx;
4938	delete [] tempV;
4939	return;
4940	}
4941
4942	void
4943	G4HEInelastic::QuasiElasticScattering(G4bool &successful,
4944	G4HEVector pv[],
4945	G4int &vecLen,
4946	G4double &excitationEnergyGNP,
4947	G4double &excitationEnergyDTA,
4948	G4HEVector incidentParticle,
4949	G4HEVector targetParticle,
4950	G4double atomicWeight,
4951	G4double atomicNumber )
4952	{
4953	// if the Cascading or Resonance - model fails, we try this,
4954	// QuasiElasticScattering.
4955	//
4956	// All quantities on the G4HEVector Array pv are in GeV- units.
4957
4958	G4int protonCode = Proton.getCode();
4959	G4String mesonType = PionPlus.getType();
4960	G4String baryonType = Proton.getType();
4961	G4String antiBaryonType= AntiProton.getType();
4962
4963	G4double targetMass = targetParticle.getMass();
4964
4965	G4double incidentMass = incidentParticle.getMass();
4966	G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
4967	G4double incidentEnergy = incidentParticle.getEnergy();
4968	G4double incidentKineticEnergy = incidentParticle.getKineticEnergy();
4969	G4String incidentType = incidentParticle.getType();
4970	// G4double incidentTOF = incidentParticle.getTOF();
4971	G4double incidentTOF = 0.;
4972
4973	// some local variables
4974
4975	G4int i;
4976
4977	if(verboseLevel > 1)
4978	G4cout << " G4HEInelastic::QuasiElasticScattering " << G4endl;
4979
4980	if (incidentTotalMomentum < 0.01 \|\| vecLen < 2 )
4981	{
4982	successful = false;
4983	return;
4984	}
4985	G4double centerOfMassEnergy = std::sqrt( sqr(incidentMass) + sqr(targetMass)
4986	+2.targetMassincidentEnergy);
4987
4988	G4HEVector pvI = incidentParticle; // for the incident particle
4989	pvI.setSide( 1 );
4990
4991	G4HEVector pvT = targetParticle; // for the target particle
4992	pvT.setMomentumAndUpdate( 0.0, 0.0, 0.0 );
4993	pvT.setSide( -1 );
4994	pvT.setTOF( -1.);
4995
4996	G4HEVector* pvmx = new G4HEVector[3];
4997
4998	if (atomicWeight > 1.5) // for the following case better use ElasticScattering.
4999	{
5000	if ( (pvI.getCode() == pv[0].getCode() )
5001	&& (pvT.getCode() == pv[1].getCode() )
5002	&& (excitationEnergyGNP < 0.001)
5003	&& (excitationEnergyDTA < 0.001) )
5004	{
5005	successful = false;
5006	return;
5007	}
5008	}
5009
5010	pv[0].setSide( 1 );
5011	pv[0].setFlag( false );
5012	pv[0].setTOF( incidentTOF);
5013	pv[0].setMomentumAndUpdate( incidentParticle.getMomentum() );
5014	pv[1].setSide( -1 );
5015	pv[1].setFlag( false );
5016	pv[1].setTOF( incidentTOF);
5017	pv[1].setMomentumAndUpdate(targetParticle.getMomentum() );
5018
5019	if ( (incidentTotalMomentum > 0.1) && (centerOfMassEnergy > 0.01) )
5020	{
5021	if ( pv[1].getType() == mesonType )
5022	{
5023	if (G4UniformRand() < 0.5)
5024	pv[1].setDefinition( "Proton");
5025	else
5026	pv[1].setDefinition( "Neutron");
5027	}
5028	pvmx[0].Add( pvI, pvT );
5029	pvmx[1].Lor( pvI, pvmx[0] );
5030	pvmx[2].Lor( pvT, pvmx[0] );
5031	G4double pin = pvmx[1].Length();
5032	G4double bvalue = Amax(0.01 , 4.225+1.795*std::log(incidentTotalMomentum));
5033	G4double pf = sqr( sqr(centerOfMassEnergy) + sqr(pv[1].getMass()) - sqr(pv[0].getMass()))
5034	- 4 * sqr(centerOfMassEnergy) * sqr(pv[1].getMass());
5035	if ( pf < 0.001)
5036	{
5037	successful = false;
5038	return;
5039	}
5040	pf = std::sqrt(pf)/(2.*centerOfMassEnergy);
5041	G4double btrang = 4. * bvalue * pin * pf;
5042	G4double exindt = -1.;
5043	if (btrang < 46.) exindt += std::exp(-btrang);
5044	G4double tdn = std::log(1. + G4UniformRand()*exindt)/btrang;
5045	G4double ctet = Amax( -1., Amin(1., 1. + 2.*tdn));
5046	G4double stet = std::sqrt((1.-ctet)*(1.+ctet));
5047	G4double phi = twopi * G4UniformRand();
5048	pv[0].setMomentumAndUpdate( pfstetstd::sin(phi),
5049	pfstetstd::cos(phi),
5050	pf*ctet );
5051	pv[1].SmulAndUpdate( pv[0], -1.);
5052
5053	for (i = 0; i < 2; i++)
5054	{
5055	pv[i].Lor( pv[i], pvmx[4] );
5056	pv[i].Defs1( pv[i], pvI );
5057	if (atomicWeight > 1.5)
5058	{
5059	G4double ekin = pv[i].getKineticEnergy()
5060	- 0.025((atomicWeight-1.)/120.)std::exp(-(atomicWeight-1.)/120.)
5061	(1. + 0.5normal());
5062	ekin = Amax(0.0001, ekin);
5063	pv[i].setKineticEnergyAndUpdate( ekin );
5064	}
5065	}
5066	}
5067	vecLen = 2;
5068
5069	// add black track particles
5070	// the total number of particles produced is restricted to 198
5071	// this may have influence on very high energies
5072
5073	if (verboseLevel > 1)
5074	G4cout << " Evaporation " << atomicWeight << " " <<
5075	excitationEnergyGNP << " " << excitationEnergyDTA << G4endl;
5076
5077	if( atomicWeight > 1.5 )
5078	{
5079
5080	G4double sprob, cost, sint, ekin2, ran, pp, eka;
5081	G4double ekin, cfa, ekin1, phi, pvMass, pvEnergy;
5082	G4int spall(0), nbl(0);
5083	// sprob is the probability of self-absorption in heavy molecules
5084
5085	sprob = 0.;
5086	cfa = 0.025((atomicWeight-1.)/120.)std::exp(-(atomicWeight-1.)/120.);
5087	// first add protons and neutrons
5088
5089	if( excitationEnergyGNP >= 0.001 )
5090	{
5091	// nbl = number of proton/neutron black track particles
5092	// tex is their total kinetic energy (GeV)
5093
5094	nbl = Poisson( excitationEnergyGNP/0.02);
5095	if( nbl > atomicWeight ) nbl = (int)(atomicWeight);
5096	if (verboseLevel > 1)
5097	G4cout << " evaporation " << nbl << " " << sprob << G4endl;
5098	spall = 0;
5099	if( nbl > 0)
5100	{
5101	ekin = excitationEnergyGNP/nbl;
5102	ekin2 = 0.0;
5103	for( i=0; i<nbl; i++ )
5104	{
5105	if( G4UniformRand() < sprob ) continue;
5106	if( ekin2 > excitationEnergyGNP) break;
5107	ran = G4UniformRand();
5108	ekin1 = -ekinstd::log(ran) - cfa(1.0+0.5*normal());
5109	if (ekin1 < 0) ekin1 = -0.010*std::log(ran);
5110	ekin2 += ekin1;
5111	if( ekin2 > excitationEnergyGNP)
5112	ekin1 = Amax( 1.0e-6, excitationEnergyGNP-(ekin2-ekin1) );
5113	if( G4UniformRand() > (1.0-atomicNumber/(atomicWeight)))
5114	pv[vecLen].setDefinition( "Proton");
5115	else
5116	pv[vecLen].setDefinition( "Neutron");
5117	spall++;
5118	cost = G4UniformRand() * 2.0 - 1.0;
5119	sint = std::sqrt(std::fabs(1.0-cost*cost));
5120	phi = twopi * G4UniformRand();
5121	pv[vecLen].setFlag( true ); // true is the same as IPA(i)<0
5122	pv[vecLen].setSide( -4 );
5123	pvMass = pv[vecLen].getMass();
5124	pv[vecLen].setTOF( 1.0 );
5125	pvEnergy = ekin1 + pvMass;
5126	pp = std::sqrt( std::fabs( pvEnergypvEnergy - pvMasspvMass ) );
5127	pv[vecLen].setMomentumAndUpdate( ppsintstd::sin(phi),
5128	ppsintstd::cos(phi),
5129	pp*cost );
5130	if (verboseLevel > 1) pv[vecLen].Print(vecLen);
5131	vecLen++;
5132	}
5133	if( (atomicWeight >= 10.0 ) && (incidentKineticEnergy <= 2.0) )
5134	{
5135	G4int ika, kk = 0;
5136	eka = incidentKineticEnergy;
5137	if( eka > 1.0 )eka *= eka;
5138	eka = Amax( 0.1, eka );
5139	ika = G4int(3.6std::exp((atomicNumberatomicNumber
5140	/atomicWeight-35.56)/6.45)/eka);
5141	if( ika > 0 )
5142	{
5143	for( i=(vecLen-1); i>=0; i-- )
5144	{
5145	if( (pv[i].getCode() == protonCode) && pv[i].getFlag() )
5146	{
5147	pv[i].setDefinition( "Neutron" );
5148	if (verboseLevel > 1) pv[i].Print(i);
5149	if( ++kk > ika ) break;
5150	}
5151	}
5152	}
5153	}
5154	}
5155	}
5156
5157	// finished adding proton/neutron black track particles
5158	// now, try to add deuterons, tritons and alphas
5159
5160	if( excitationEnergyDTA >= 0.001 )
5161	{
5162	nbl = (G4int)(2.*std::log(atomicWeight));
5163
5164	// nbl is the number of deutrons, tritons, and alphas produced
5165
5166	if( nbl > 0 )
5167	{
5168	ekin = excitationEnergyDTA/nbl;
5169	ekin2 = 0.0;
5170	for( i=0; i<nbl; i++ )
5171	{
5172	if( G4UniformRand() < sprob ) continue;
5173	if( ekin2 > excitationEnergyDTA) break;
5174	ran = G4UniformRand();
5175	ekin1 = -ekinstd::log(ran)-cfa(1.+0.5*normal());
5176	if( ekin1 < 0.0 ) ekin1 = -0.010*std::log(ran);
5177	ekin2 += ekin1;
5178	if( ekin2 > excitationEnergyDTA)
5179	ekin1 = Amax( 1.0e-6, excitationEnergyDTA-(ekin2-ekin1));
5180	cost = G4UniformRand()*2.0 - 1.0;
5181	sint = std::sqrt(std::fabs(1.0-cost*cost));
5182	phi = twopi*G4UniformRand();
5183	ran = G4UniformRand();
5184	if( ran <= 0.60 )
5185	pv[vecLen].setDefinition( "Deuteron");
5186	else if (ran <= 0.90)
5187	pv[vecLen].setDefinition( "Triton");
5188	else
5189	pv[vecLen].setDefinition( "Alpha");
5190	spall += (int)(pv[vecLen].getMass() * 1.066);
5191	if( spall > atomicWeight ) break;
5192	pv[vecLen].setFlag( true ); // true is the same as IPA(i)<0
5193	pv[vecLen].setSide( -4 );
5194	pvMass = pv[vecLen].getMass();
5195	pv[vecLen].setTOF( 1.0 );
5196	pvEnergy = pvMass + ekin1;
5197	pp = std::sqrt( std::fabs( pvEnergypvEnergy - pvMasspvMass ) );
5198	pv[vecLen].setMomentumAndUpdate( ppsintstd::sin(phi),
5199	ppsintstd::cos(phi),
5200	pp*cost );
5201	if (verboseLevel > 1) pv[vecLen].Print(vecLen);
5202	vecLen++;
5203	}
5204	}
5205	}
5206	}
5207
5208	// Calculate time delay for nuclear reactions
5209
5210	G4double tof = incidentTOF;
5211	if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0)
5212	&& (incidentKineticEnergy <= 0.2) )
5213	tof -= 500.0 * std::exp(-incidentKineticEnergy /0.04) * std::log( G4UniformRand() );
5214	for ( i=0; i < vecLen; i++)
5215	{
5216
5217	pv[i].setTOF ( tof );
5218	// vec[i].SetTOF ( tof );
5219	}
5220
5221	for(i=0; i<vecLen; i++)
5222	{
5223	if(pv[i].getName() == "KaonZero" \|\| pv[i].getName() == "AntiKaonZero")
5224	{
5225	pvmx[0] = pv[i];
5226	if(G4UniformRand() < 0.5) pv[i].setDefinition("KaonZeroShort");
5227	else pv[i].setDefinition("KaonZeroLong");
5228	pv[i].setMomentumAndUpdate(pvmx[0].getMomentum());
5229	}
5230	}
5231
5232	successful = true;
5233	delete [] pvmx;
5234	return;
5235	}
5236
5237	void
5238	G4HEInelastic::ElasticScattering(G4bool &successful,
5239	G4HEVector pv[],
5240	G4int &vecLen,
5241	G4HEVector incidentParticle,
5242	G4double atomicWeight,
5243	G4double /* atomicNumber*/)
5244	{
5245	if(verboseLevel > 1)
5246	G4cout << " G4HEInelastic::ElasticScattering " << G4endl;
5247
5248	G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
5249	if (verboseLevel > 1)
5250	G4cout << "DoIt: Incident particle momentum="
5251	<< incidentTotalMomentum << " GeV" << G4endl;
5252	if (incidentTotalMomentum < 0.01)
5253	{
5254	successful = false;
5255	return;
5256	}
5257	if (atomicWeight < 0.5)
5258	{
5259	successful = false;
5260	return;
5261	}
5262	pv[0] = incidentParticle;
5263	vecLen = 1;
5264
5265	G4double aa, bb, cc, dd, rr;
5266	if (atomicWeight <= 62.)
5267	{
5268	aa = std::pow(atomicWeight, 1.63);
5269	bb = 14.5*std::pow(atomicWeight, 0.66);
5270	cc = 1.4*std::pow(atomicWeight, 0.33);
5271	dd = 10.;
5272	}
5273	else
5274	{
5275	aa = std::pow(atomicWeight, 1.33);
5276	bb = 60.*std::pow(atomicWeight, 0.33);
5277	cc = 0.4*std::pow(atomicWeight, 0.40);
5278	dd = 10.;
5279	}
5280	aa = aa/bb;
5281	cc = cc/dd;
5282	G4double ran = G4UniformRand();
5283	rr = (aa + cc)*ran;
5284	if (verboseLevel > 1)
5285	{
5286	G4cout << "ElasticScattering: aa,bb,cc,dd,rr" << G4endl;
5287	G4cout << aa << " " << bb << " " << cc << " " << dd << " "
5288	<< rr << G4endl;
5289	}
5290	G4double t1 = -std::log(ran)/bb;
5291	G4double t2 = -std::log(ran)/dd;
5292	if (verboseLevel > 1) {
5293	G4cout << "t1,fctcos " << t1 << " " << fctcos(t1, aa, bb, cc, dd, rr)
5294	<< G4endl;
5295	G4cout << "t2,fctcos " << t2 << " " << fctcos(t2, aa, bb, cc, dd, rr)
5296	<< G4endl;
5297	}
5298	G4double eps = 0.001;
5299	G4int ind1 = 10;
5300	G4double t;
5301	G4int ier1;
5302	ier1 = rtmi(&t, t1, t2, eps, ind1, aa, bb, cc, dd, rr);
5303	if (verboseLevel > 1) {
5304	G4cout << "From rtmi, ier1=" << ier1 << G4endl;
5305	G4cout << "t, fctcos " << t << " " << fctcos(t, aa, bb, cc, dd, rr)
5306	<< G4endl;
5307	}
5308	if (ier1 != 0) t = 0.25(3.t1 + t2);
5309	if (verboseLevel > 1)
5310	G4cout << "t, fctcos " << t << " " << fctcos(t, aa, bb, cc, dd, rr)
5311	<< G4endl;
5312
5313	G4double phi = G4UniformRand()*twopi;
5314	rr = 0.5*t/sqr(incidentTotalMomentum);
5315	if (rr > 1.) rr = 0.;
5316	if (verboseLevel > 1)
5317	G4cout << "rr=" << rr << G4endl;
5318	G4double cost = 1. - rr;
5319	G4double sint = std::sqrt(Amax(rr*(2. - rr), 0.));
5320	if (verboseLevel > 1)
5321	G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
5322	// Scattered particle referred to axis of incident particle
5323	G4HEVector pv0;
5324	G4HEVector pvI;
5325	pvI.setMass( incidentParticle.getMass() );
5326	pvI.setMomentum( incidentParticle.getMomentum() );
5327	pvI.SmulAndUpdate( pvI, 1. );
5328	pv0.setMass( pvI.getMass() );
5329
5330	pv0.setMomentumAndUpdate( incidentTotalMomentum * sint * std::sin(phi),
5331	incidentTotalMomentum * sint * std::cos(phi),
5332	incidentTotalMomentum * cost );
5333	pv0.Defs1( pv0, pvI );
5334
5335	successful = true;
5336	return;
5337	}
5338
5339
5340	G4int
5341	G4HEInelastic::rtmi(G4double *x, G4double xli, G4double xri, G4double eps,
5342	G4int iend,
5343	G4double aa, G4double bb, G4double cc, G4double dd,
5344	G4double rr)
5345	{
5346	G4int ier = 0;
5347	G4double xl = xli;
5348	G4double xr = xri;
5349	*x = xl;
5350	G4double tol = *x;
5351	G4double f = fctcos(tol, aa, bb, cc, dd, rr);
5352	if (f == 0.) return ier;
5353	G4double fl, fr;
5354	fl = f;
5355	*x = xr;
5356	tol = *x;
5357	f = fctcos(tol, aa, bb, cc, dd, rr);
5358	if (f == 0.) return ier;
5359	fr = f;
5360
5361	// Error return in case of wrong input data
5362	if (fl*fr >= 0.)
5363	{
5364	ier = 2;
5365	return ier;
5366	}
5367
5368	// Basic assumption fl*fr less than 0 is satisfied.
5369	// Generate tolerance for function values.
5370	G4int i = 0;
5371	G4double tolf = 100.*eps;
5372
5373	// Start iteration loop
5374
5375	label4: // <-------------
5376	i++;
5377
5378	// Start bisection loop
5379
5380	for (G4int k = 1; k <= iend; k++)
5381	{
5382	x = 0.5(xl + xr);
5383	tol = *x;
5384	f = fctcos(tol, aa, bb, cc, dd, rr);
5385	if (f == 0.) return 0;
5386	if (f*fr < 0.)
5387	{ // Interchange xl and xr in order to get the
5388	tol = xl; // same sign in f and fr
5389	xl = xr;
5390	xr = tol;
5391	tol = fl;
5392	fl = fr;
5393	fr = tol;
5394	}
5395	tol = f - fl;
5396	G4double a = f*tol;
5397	a = a + a;
5398	if (a < fr*(fr - fl) && i <= iend) goto label17;
5399	xr = *x;
5400	fr = f;
5401
5402	// Test on satisfactory accuracy in bisection loop
5403	tol = eps;
5404	a = std::fabs(xr);
5405	if (a > 1.) tol = tol*a;
5406	if (std::fabs(xr - xl) <= tol && std::fabs(fr - fl) <= tolf) goto label14;
5407	}
5408	// End of bisection loop
5409
5410	// No convergence after iend iteration steps followed by iend
5411	// successive steps of bisection or steadily increasing function
5412	// values at right bounds. Error return.
5413	ier = 1;
5414
5415	label14: // <---------------
5416	if (std::fabs(fr) > std::fabs(fl))
5417	{
5418	*x = xl;
5419	f = fl;
5420	}
5421	return ier;
5422
5423	// Computation of iterated x-value by inverse parabolic interp
5424	label17: // <---------------
5425	G4double a = fr - f;
5426	G4double dx = (x - xl)fl(1. + f(a - tol)/(a*(fr - fl)))/tol;
5427	G4double xm = *x;
5428	G4double fm = f;
5429	*x = xl - dx;
5430	tol = *x;
5431	f = fctcos(tol, aa, bb, cc, dd, rr);
5432	if (f == 0.) return ier;
5433
5434	// Test on satisfactory accuracy in iteration loop
5435	tol = eps;
5436	a = std::fabs(*x);
5437	if (a > 1) tol = tol*a;
5438	if (std::fabs(dx) <= tol && std::fabs(f) <= tolf) return ier;
5439
5440	// Preparation of next bisection loop
5441	if (f*fl < 0.)
5442	{
5443	xr = *x;
5444	fr = f;
5445	}
5446	else
5447	{
5448	xl = *x;
5449	fl = f;
5450	xr = xm;
5451	fr = fm;
5452	}
5453	goto label4;
5454	}
5455
5456
5457	// Test function for root-finder
5458
5459	G4double
5460	G4HEInelastic::fctcos(G4double t, G4double aa, G4double bb, G4double cc,
5461	G4double dd, G4double rr)
5462	{
5463	const G4double expxl = -82.;
5464	const G4double expxu = 82.;
5465
5466	G4double test1 = -bb*t;
5467	if (test1 > expxu) test1 = expxu;
5468	if (test1 < expxl) test1 = expxl;
5469
5470	G4double test2 = -dd*t;
5471	if (test2 > expxu) test2 = expxu;
5472	if (test2 < expxl) test2 = expxl;
5473
5474	return aastd::exp(test1) + ccstd::exp(test2) - rr;
5475	}
5476
5477	G4double G4HEInelastic::NBodyPhaseSpace
5478	( const G4double totalEnergy, // MeV
5479	const G4bool constantCrossSection,
5480	G4HEVector vec[],
5481	G4int& vecLen )
5482	{
5483	// derived from original FORTRAN code PHASP by H. Fesefeldt (02-Dec-1986)
5484	// Returns the weight of the event
5485
5486	G4int i;
5487
5488	const G4double expxu = std::log(FLT_MAX); // upper bound for arg. of exp
5489	const G4double expxl = -expxu; // lower bound for arg. of exp
5490
5491	if( vecLen < 2 ) {
5492	G4cerr << "*** Error in G4HEInelastic::GenerateNBodyEvent" << G4endl;
5493	G4cerr << " number of particles < 2" << G4endl;
5494	G4cerr << "totalEnergy = " << totalEnergy << ", vecLen = "
5495	<< vecLen << G4endl;
5496	return -1.0;
5497	}
5498
5499	G4double* mass = new G4double [vecLen]; // mass of each particle
5500	G4double* energy = new G4double [vecLen]; // total energy of each particle
5501	G4double** pcm; // pcm is an array with 3 rows and vecLen columns
5502	pcm = new G4double* [3];
5503	for( i=0; i<3; ++i )pcm[i] = new G4double [vecLen];
5504
5505	G4double totalMass = 0.0;
5506	G4double* sm = new G4double [vecLen];
5507
5508	for( i=0; i<vecLen; ++i ) {
5509	mass[i] = vec[i].getMass();
5510	vec[i].setMomentum( 0.0, 0.0, 0.0 );
5511	pcm[0][i] = 0.0; // x-momentum of i-th particle
5512	pcm[1][i] = 0.0; // y-momentum of i-th particle
5513	pcm[2][i] = 0.0; // z-momentum of i-th particle
5514	energy[i] = mass[i]; // total energy of i-th particle
5515	totalMass += mass[i];
5516	sm[i] = totalMass;
5517	}
5518
5519	if( totalMass >= totalEnergy ) {
5520	if (verboseLevel > 1) {
5521	G4cout << "*** Error in G4HEInelastic::GenerateNBodyEvent" << G4endl;
5522	G4cout << " total mass (" << totalMass << ") >= total energy ("
5523	<< totalEnergy << ")" << G4endl;
5524	}
5525	delete [] mass;
5526	delete [] energy;
5527	for( i=0; i<3; ++i )delete [] pcm[i];
5528	delete [] pcm;
5529	delete [] sm;
5530	return -1.0;
5531	}
5532
5533	G4double kineticEnergy = totalEnergy - totalMass;
5534	G4double* emm = new G4double [vecLen];
5535	emm[0] = mass[0];
5536	if( vecLen > 3 ) { // the random numbers are sorted
5537	G4double* ran = new G4double [vecLen];
5538	for( i=0; i<vecLen; ++i )ran[i] = G4UniformRand();
5539	for( i=0; i<vecLen-1; ++i ) {
5540	for( G4int j=vecLen-1; j > i; --j ) {
5541	if( ran[i] > ran[j] ) {
5542	G4double temp = ran[i];
5543	ran[i] = ran[j];
5544	ran[j] = temp;
5545	}
5546	}
5547	}
5548	for( i=1; i<vecLen; ++i )emm[i] = ran[i-1]*kineticEnergy + sm[i];
5549	delete [] ran;
5550	} else {
5551	emm[1] = G4UniformRand()*kineticEnergy + sm[1];
5552	}
5553	emm[vecLen-1] = totalEnergy;
5554
5555	// Weight is the sum of logarithms of terms instead of the product of terms
5556
5557	G4bool lzero = true;
5558	G4double wtmax = 0.0;
5559	if( constantCrossSection ) { // this is KGENEV=1 in PHASP
5560	G4double emmax = kineticEnergy + mass[0];
5561	G4double emmin = 0.0;
5562	for( i=1; i<vecLen; ++i ) {
5563	emmin += mass[i-1];
5564	emmax += mass[i];
5565	G4double wtfc = 0.0;
5566	if( emmax*emmax > 0.0 ) {
5567	G4double arg = emmax*emmax
5568	+ (emminemmin-mass[i]mass[i])(emminemmin-mass[i]mass[i])/(emmaxemmax)
5569	- 2.0(emminemmin+mass[i]*mass[i]);
5570	if( arg > 0.0 )wtfc = 0.5*std::sqrt( arg );
5571	}
5572	if( wtfc == 0.0 ) {
5573	lzero = false;
5574	break;
5575	}
5576	wtmax += std::log( wtfc );
5577	}
5578	if( lzero )
5579	wtmax = -wtmax;
5580	else
5581	wtmax = expxu;
5582	} else {
5583	wtmax = std::log( std::pow( kineticEnergy, vecLen-2 ) *
5584	pi * std::pow( twopi, vecLen-2 ) / Factorial(vecLen-2) );
5585	}
5586	lzero = true;
5587	G4double* pd = new G4double [vecLen-1];
5588	for( i=0; i<vecLen-1; ++i ) {
5589	pd[i] = 0.0;
5590	if( emm[i+1]*emm[i+1] > 0.0 ) {
5591	G4double arg = emm[i+1]*emm[i+1]
5592	+ (emm[i]emm[i]-mass[i+1]mass[i+1])(emm[i]emm[i]-mass[i+1]*mass[i+1])
5593	/(emm[i+1]*emm[i+1])
5594	- 2.0(emm[i]emm[i]+mass[i+1]*mass[i+1]);
5595	if( arg > 0.0 )pd[i] = 0.5*std::sqrt( arg );
5596	}
5597	if( pd[i] == 0.0 )
5598	lzero = false;
5599	else
5600	wtmax += std::log( pd[i] );
5601	}
5602	G4double weight = 0.0; // weight is returned by GenerateNBodyEvent
5603	if( lzero )weight = std::exp( Amax(Amin(wtmax,expxu),expxl) );
5604
5605	G4double bang, cb, sb, s0, s1, s2, c, s, esys, a, b, gama, beta;
5606	pcm[0][0] = 0.0;
5607	pcm[1][0] = pd[0];
5608	pcm[2][0] = 0.0;
5609	for( i=1; i<vecLen; ++i ) {
5610	pcm[0][i] = 0.0;
5611	pcm[1][i] = -pd[i-1];
5612	pcm[2][i] = 0.0;
5613	bang = twopi*G4UniformRand();
5614	cb = std::cos(bang);
5615	sb = std::sin(bang);
5616	c = 2.0*G4UniformRand() - 1.0;
5617	s = std::sqrt( std::fabs( 1.0-c*c ) );
5618	if( i < vecLen-1 ) {
5619	esys = std::sqrt(pd[i]pd[i] + emm[i]emm[i]);
5620	beta = pd[i]/esys;
5621	gama = esys/emm[i];
5622	for( G4int j=0; j<=i; ++j ) {
5623	s0 = pcm[0][j];
5624	s1 = pcm[1][j];
5625	s2 = pcm[2][j];
5626	energy[j] = std::sqrt( s0s0 + s1s1 + s2s2 + mass[j]mass[j] );
5627	a = s0c - s1s; // rotation
5628	pcm[1][j] = s0s + s1c;
5629	b = pcm[2][j];
5630	pcm[0][j] = acb - bsb;
5631	pcm[2][j] = asb + bcb;
5632	pcm[1][j] = gama(pcm[1][j] + betaenergy[j]);
5633	}
5634	} else {
5635	for( G4int j=0; j<=i; ++j ) {
5636	s0 = pcm[0][j];
5637	s1 = pcm[1][j];
5638	s2 = pcm[2][j];
5639	energy[j] = std::sqrt( s0s0 + s1s1 + s2s2 + mass[j]mass[j] );
5640	a = s0c - s1s; // rotation
5641	pcm[1][j] = s0s + s1c;
5642	b = pcm[2][j];
5643	pcm[0][j] = acb - bsb;
5644	pcm[2][j] = asb + bcb;
5645	}
5646	}
5647	}
5648	G4double pModule;
5649	for( i=0; i<vecLen; ++i ) {
5650	kineticEnergy = energy[i] - mass[i];
5651	pModule = std::sqrt( sqr(kineticEnergy) + 2kineticEnergymass[i] );
5652	vec[i].setMomentum( pcm[0][i]/pModule,
5653	pcm[1][i]/pModule,
5654	pcm[2][i]/pModule );
5655	vec[i].setKineticEnergyAndUpdate( kineticEnergy );
5656	}
5657	delete [] mass;
5658	delete [] energy;
5659	for( i=0; i<3; ++i )delete [] pcm[i];
5660	delete [] pcm;
5661	delete [] emm;
5662	delete [] sm;
5663	delete [] pd;
5664	return weight;
5665	}
5666
5667	G4double
5668	G4HEInelastic::gpdk( G4double a, G4double b, G4double c )
5669	{
5670	if( a == 0.0 )
5671	{
5672	return 0.0;
5673	}
5674	else
5675	{
5676	G4double arg = aa+(bb-cc)(bb-cc)/(aa)-2.0(bb+cc);
5677	if( arg <= 0.0 )
5678	{
5679	return 0.0;
5680	}
5681	else
5682	{
5683	return 0.5*std::sqrt(std::fabs(arg));
5684	}
5685	}
5686	}
5687
5688
5689	G4double
5690	G4HEInelastic::NBodyPhaseSpace(G4int npart, G4HEVector pv[],
5691	G4double wmax, G4double wfcn,
5692	G4int maxtrial, G4int ntrial)
5693	{ ntrial = 0;
5694	G4double wps(0);
5695	while ( ntrial < maxtrial)
5696	{ ntrial++;
5697	G4int i, j;
5698	G4int nrn = 3*(npart-2)-4;
5699	G4double *ranarr = new G4double[nrn];
5700	for (i=0;i<nrn;i++) ranarr[i]=G4UniformRand();
5701	G4int nrnp = npart-4;
5702	if(nrnp > 1) QuickSort( ranarr, 0 , nrnp-1 );
5703	G4HEVector pvcms;
5704	pvcms.Add(pv[0],pv[1]);
5705	pvcms.Smul( pvcms, -1.);
5706	G4double rm = 0.;
5707	for (i=2;i<npart;i++) rm += pv[i].getMass();
5708	G4double rm1 = pvcms.getMass() - rm;
5709	rm -= pv[2].getMass();
5710	wps = (npart-3)std::pow(rm1/sqr(twopi), npart-4)/(4pi*pvcms.getMass());
5711	for (i=3; (i=npart-1);i++) wps /= i-2; // @@@@@@@@@@ bug @@@@@@@@@
5712	G4double xxx = rm1/sqr(twopi);
5713	for (i=1; (i=npart-4); i++) wps /= xxx/i; // @@@@@@@@@@ bug @@@@@@@@@
5714	wps /= (4pipvcms.getMass());
5715	G4double p2,cost,sint,phi;
5716	j = 1;
5717	while (j)
5718	{ j++;
5719	rm -= pv[j+1].getMass();
5720	if(j == npart-2) break;
5721	G4double rmass = rm + rm1*ranarr[npart-j-1];
5722	p2 = Alam(sqr(pvcms.getMass()), sqr(pv[j].getMass()),
5723	sqr(rmass))/(4.*sqr(pvcms.getMass()));
5724	cost = 1. - 2.ranarr[npart+2j-9];
5725	sint = std::sqrt(1.-cost*cost);
5726	phi = twopiranarr[npart+2j-8];
5727	p2 = std::sqrt( Amax(0., p2));
5728	wps *= p2;
5729	pv[j].setMomentumAndUpdate( p2sintstd::sin(phi), p2sintstd::cos(phi),p2*cost);
5730	pv[j].Lor(pv[j], pvcms);
5731	pvcms.Add3( pvcms, pv[j] );
5732	pvcms.setEnergy(pvcms.getEnergy()-pv[j].getEnergy());
5733	pvcms.setMass( std::sqrt(sqr(pvcms.getEnergy()) - sqr(pvcms.Length())));
5734	}
5735	p2 = Alam(sqr(pvcms.getMass()), sqr(pv[j].getMass()),
5736	sqr(rm))/(4.*sqr(pvcms.getMass()));
5737	cost = 1. - 2.ranarr[npart+2j-9];
5738	sint = std::sqrt(1.-cost*cost);
5739	phi = twopiranarr[npart+2j-8];
5740	p2 = std::sqrt( Amax(0. , p2));
5741	wps *= p2;
5742	pv[j].setMomentumAndUpdate( p2sintstd::sin(phi), p2sintstd::cos(phi), p2*cost);
5743	pv[j+1].setMomentumAndUpdate( -p2sintstd::sin(phi), -p2sintstd::cos(phi), -p2*cost);
5744	pv[j].Lor( pv[j], pvcms );
5745	pv[j+1].Lor( pv[j+1], pvcms );
5746	wfcn = CalculatePhaseSpaceWeight( npart );
5747	G4double wt = wps * wfcn;
5748	if (wt > wmax)
5749	{ wmax = wt;
5750	G4cout << "maximum weight changed to " << wmax << G4endl;
5751	}
5752	wt = wt/wmax;
5753	if (G4UniformRand() < wt) break;
5754	}
5755	return wps;
5756	}
5757
5758
5759	void
5760	G4HEInelastic::QuickSort(G4double arr[], const G4int lidx, const G4int ridx)
5761	{ // sorts the Array arr[] in ascending order
5762	G4double buffer;
5763	G4int k, e, mid;
5764	if(lidx>=ridx) return;
5765	mid = (int)((lidx+ridx)/2.);
5766	buffer = arr[lidx];
5767	arr[lidx]= arr[mid];
5768	arr[mid] = buffer;
5769	e = lidx;
5770	for (k=lidx+1;k<=ridx;k++)
5771	if (arr[k] < arr[lidx])
5772	{ e++;
5773	buffer = arr[e];
5774	arr[e] = arr[k];
5775	arr[k] = buffer;
5776	}
5777	buffer = arr[lidx];
5778	arr[lidx]= arr[e];
5779	arr[e] = buffer;
5780	QuickSort(arr, lidx, e-1);
5781	QuickSort(arr, e+1 , ridx);
5782	return;
5783	}
5784
5785	G4double
5786	G4HEInelastic::Alam( G4double a, G4double b, G4double c)
5787	{ return aa + bb + cc - 2.ab - 2.ac -2.b*c;
5788	}
5789
5790	G4double
5791	G4HEInelastic::CalculatePhaseSpaceWeight( G4int /* npart */)
5792	{ G4double wfcn = 1.;
5793	return wfcn;
5794	}
5795
5796
5797
5798
5799
5800
5801

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

Download in other formats: