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

Last change on this file since 893 was 819, checked in by garnier, 17 years ago

import all except CVS

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