Context Navigation

G4HEInelastic.cc @ 1347

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

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

Context Navigation

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

Download in other formats: