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