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