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

Last change on this file since 1350 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4

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