Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGReaction.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 // $Id$
27 //
28 
29 #include <iostream>
30 
31 #include "G4RPGReaction.hh"
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "Randomize.hh"
35 
36 G4bool G4RPGReaction::
37 ReactionStage(const G4HadProjectile* /*originalIncident*/,
38  G4ReactionProduct& /*modifiedOriginal*/,
39  G4bool& /*incidentHasChanged*/,
40  const G4DynamicParticle* /*originalTarget*/,
41  G4ReactionProduct& /*targetParticle*/,
42  G4bool& /*targetHasChanged*/,
43  const G4Nucleus& /*targetNucleus*/,
44  G4ReactionProduct& /*currentParticle*/,
45  G4FastVector<G4ReactionProduct,256>& /*vec*/,
46  G4int& /*vecLen*/,
47  G4bool /*leadFlag*/,
48  G4ReactionProduct& /*leadingStrangeParticle*/)
49 {
50  G4cout << " G4RPGReactionStage must be overridden in a derived class "
51  << G4endl;
52  return false;
53 }
54 
55 
56 void G4RPGReaction::
57 AddBlackTrackParticles(const G4double epnb, // GeV
58  const G4int npnb,
59  const G4double edta, // GeV
60  const G4int ndta,
61  const G4ReactionProduct& modifiedOriginal,
62  G4int PinNucleus,
63  G4int NinNucleus,
64  const G4Nucleus& targetNucleus,
65  G4FastVector<G4ReactionProduct,256>& vec,
66  G4int& vecLen)
67 {
68  // derived from original FORTRAN code in GENXPT and TWOCLU by H. Fesefeldt
69  //
70  // npnb is number of proton/neutron black track particles
71  // ndta is the number of deuterons, tritons, and alphas produced
72  // epnb is the kinetic energy available for proton/neutron black track particles
73  // edta is the kinetic energy available for deuteron/triton/alpha particles
74 
75  G4ParticleDefinition* aProton = G4Proton::Proton();
76  G4ParticleDefinition* aNeutron = G4Neutron::Neutron();
77  G4ParticleDefinition* aDeuteron = G4Deuteron::Deuteron();
78  G4ParticleDefinition* aTriton = G4Triton::Triton();
79  G4ParticleDefinition* anAlpha = G4Alpha::Alpha();
80 
81  const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/MeV;
82  const G4double atomicWeight = targetNucleus.GetA_asInt();
83  const G4double atomicNumber = targetNucleus.GetZ_asInt();
84 
85  const G4double ika1 = 3.6;
86  const G4double ika2 = 35.56;
87  const G4double ika3 = 6.45;
88 
89  G4int i;
90  G4double pp;
91  G4double kinetic = 0;
92  G4double kinCreated = 0;
93  // G4double cfa = 0.025*((atomicWeight-1.0)/120.0) * std::exp(-(atomicWeight-1.0)/120.0);
94  G4double remainingE = 0;
95 
96  // First add protons and neutrons to final state
97  if (npnb > 0) {
98  // G4double backwardKinetic = 0.0;
99  G4int local_npnb = npnb;
100  // DHW: does not conserve energy for (i = 0; i < npnb; ++i) if (G4UniformRand() < sprob) local_npnb--;
101  local_npnb = std::min(PinNucleus + NinNucleus , local_npnb);
102  G4double local_epnb = epnb;
103  if (ndta == 0) local_epnb += edta; // Retrieve unused kinetic energy
104  // G4double ekin = local_epnb/std::max(1,local_npnb);
105 
106  remainingE = local_epnb;
107  for (i = 0; i < local_npnb; ++i)
108  {
109  G4ReactionProduct* p1 = new G4ReactionProduct();
110  // if( backwardKinetic > local_epnb ) {
111  // delete p1;
112  // break;
113  // }
114 
115  // G4double ran = G4UniformRand();
116  // G4double kinetic = -ekin*std::log(ran) - cfa*(1.0+0.5*normal());
117  // if( kinetic < 0.0 )kinetic = -0.010*std::log(ran);
118  // backwardKinetic += kinetic;
119  // if( backwardKinetic > local_epnb )
120  // kinetic = std::max( kineticMinimum, local_epnb-(backwardKinetic-kinetic) );
121 
122  if (G4UniformRand() > (1.0-atomicNumber/atomicWeight)) {
123 
124  // Boil off a proton if there are any left, otherwise a neutron
125 
126  if (PinNucleus > 0) {
127  p1->SetDefinition( aProton );
128  PinNucleus--;
129  } else {
130  p1->SetDefinition( aNeutron );
131  NinNucleus--;
132  // } else {
133  // delete p1;
134  // break; // no nucleons left in nucleus
135  }
136  } else {
137 
138  // Boil off a neutron if there are any left, otherwise a proton
139 
140  if (NinNucleus > 0) {
141  p1->SetDefinition( aNeutron );
142  NinNucleus--;
143  } else {
144  p1->SetDefinition( aProton );
145  PinNucleus--;
146  // } else {
147  // delete p1;
148  // break; // no nucleons left in nucleus
149  }
150  }
151 
152  if (i < local_npnb - 1) {
153  kinetic = remainingE*G4UniformRand();
154  remainingE -= kinetic;
155  } else {
156  kinetic = remainingE;
157  }
158 
159  vec.SetElement( vecLen, p1 );
160  G4double cost = G4UniformRand() * 2.0 - 1.0;
161  G4double sint = std::sqrt(std::fabs(1.0-cost*cost));
162  G4double phi = twopi * G4UniformRand();
163  vec[vecLen]->SetNewlyAdded( true );
164  vec[vecLen]->SetKineticEnergy( kinetic*GeV );
165  kinCreated+=kinetic;
166  pp = vec[vecLen]->GetTotalMomentum();
167  vec[vecLen]->SetMomentum(pp*sint*std::sin(phi),
168  pp*sint*std::cos(phi),
169  pp*cost );
170  vecLen++;
171  }
172 
173  if (NinNucleus > 0) {
174  if( (atomicWeight >= 10.0) && (ekOriginal <= 2.0*GeV) )
175  {
176  G4double ekw = ekOriginal/GeV;
177  G4int ika, kk = 0;
178  if( ekw > 1.0 )ekw *= ekw;
179  ekw = std::max( 0.1, ekw );
180  ika = G4int(ika1*std::exp((atomicNumber*atomicNumber/
181  atomicWeight-ika2)/ika3)/ekw);
182  if( ika > 0 )
183  {
184  for( i=(vecLen-1); i>=0; --i )
185  {
186  if( (vec[i]->GetDefinition() == aProton) && vec[i]->GetNewlyAdded() )
187  {
188  vec[i]->SetDefinitionAndUpdateE( aNeutron ); // modified 22-Oct-97
189  PinNucleus++;
190  NinNucleus--;
191  if( ++kk > ika )break;
192  }
193  }
194  }
195  }
196  } // if (NinNucleus >0)
197  } // if (npnb > 0)
198 
199  // Next try to add deuterons, tritons and alphas to final state
200 
201  G4double ran = 0;
202  if (ndta > 0) {
203  // G4double backwardKinetic = 0.0;
204  G4int local_ndta=ndta;
205  // DHW: does not conserve energy for (i = 0; i < ndta; ++i) if (G4UniformRand() < sprob) local_ndta--;
206  G4double local_edta = edta;
207  if (npnb == 0) local_edta += epnb; // Retrieve unused kinetic energy
208  // G4double ekin = local_edta/std::max(1,local_ndta);
209 
210  remainingE = local_edta;
211  for (i = 0; i < local_ndta; ++i) {
212  G4ReactionProduct* p2 = new G4ReactionProduct();
213  // if( backwardKinetic > local_edta ) {
214  // delete p2;
215  // break;
216  // }
217 
218  // G4double ran = G4UniformRand();
219  // G4double kinetic = -ekin*std::log(ran)-cfa*(1.+0.5*normal());
220  // if( kinetic < 0.0 )kinetic = kineticFactor*std::log(ran);
221  // backwardKinetic += kinetic;
222  // if( backwardKinetic > local_edta )kinetic = local_edta-(backwardKinetic-kinetic);
223  // if( kinetic < 0.0 )kinetic = kineticMinimum;
224 
225  ran = G4UniformRand();
226  if (ran < 0.60) {
227  if (PinNucleus > 0 && NinNucleus > 0) {
228  p2->SetDefinition( aDeuteron );
229  PinNucleus--;
230  NinNucleus--;
231  } else if (NinNucleus > 0) {
232  p2->SetDefinition( aNeutron );
233  NinNucleus--;
234  } else if (PinNucleus > 0) {
235  p2->SetDefinition( aProton );
236  PinNucleus--;
237  } else {
238  delete p2;
239  break;
240  }
241  } else if (ran < 0.90) {
242  if (PinNucleus > 0 && NinNucleus > 1) {
243  p2->SetDefinition( aTriton );
244  PinNucleus--;
245  NinNucleus -= 2;
246  } else if (PinNucleus > 0 && NinNucleus > 0) {
247  p2->SetDefinition( aDeuteron );
248  PinNucleus--;
249  NinNucleus--;
250  } else if (NinNucleus > 0) {
251  p2->SetDefinition( aNeutron );
252  NinNucleus--;
253  } else if (PinNucleus > 0) {
254  p2->SetDefinition( aProton );
255  PinNucleus--;
256  } else {
257  delete p2;
258  break;
259  }
260  } else {
261  if (PinNucleus > 1 && NinNucleus > 1) {
262  p2->SetDefinition( anAlpha );
263  PinNucleus -= 2;
264  NinNucleus -= 2;
265  } else if (PinNucleus > 0 && NinNucleus > 1) {
266  p2->SetDefinition( aTriton );
267  PinNucleus--;
268  NinNucleus -= 2;
269  } else if (PinNucleus > 0 && NinNucleus > 0) {
270  p2->SetDefinition( aDeuteron );
271  PinNucleus--;
272  NinNucleus--;
273  } else if (NinNucleus > 0) {
274  p2->SetDefinition( aNeutron );
275  NinNucleus--;
276  } else if (PinNucleus > 0) {
277  p2->SetDefinition( aProton );
278  PinNucleus--;
279  } else {
280  delete p2;
281  break;
282  }
283  }
284 
285  if (i < local_ndta - 1) {
286  kinetic = remainingE*G4UniformRand();
287  remainingE -= kinetic;
288  } else {
289  kinetic = remainingE;
290  }
291 
292  vec.SetElement( vecLen, p2 );
293  G4double cost = 2.0*G4UniformRand() - 1.0;
294  G4double sint = std::sqrt(std::max(0.0,(1.0-cost*cost)));
295  G4double phi = twopi*G4UniformRand();
296  vec[vecLen]->SetNewlyAdded( true );
297  vec[vecLen]->SetKineticEnergy( kinetic*GeV );
298  kinCreated+=kinetic;
299 
300  pp = vec[vecLen]->GetTotalMomentum();
301  vec[vecLen]->SetMomentum( pp*sint*std::sin(phi),
302  pp*sint*std::cos(phi),
303  pp*cost );
304  vecLen++;
305  }
306  } // if (ndta > 0)
307 }
308 
309 
310 G4double
311 G4RPGReaction::GenerateNBodyEvent(const G4double totalEnergy, // MeV
312  const G4bool constantCrossSection,
313  G4FastVector<G4ReactionProduct,256>& vec,
314  G4int &vecLen)
315 {
316  G4int i;
317  const G4double expxu = 82.; // upper bound for arg. of exp
318  const G4double expxl = -expxu; // lower bound for arg. of exp
319 
320  if (vecLen < 2) {
321  G4cerr << "*** Error in G4RPGReaction::GenerateNBodyEvent" << G4endl;
322  G4cerr << " number of particles < 2" << G4endl;
323  G4cerr << "totalEnergy = " << totalEnergy << "MeV, vecLen = " << vecLen << G4endl;
324  return -1.0;
325  }
326 
327  G4double mass[18]; // mass of each particle
328  G4double energy[18]; // total energy of each particle
329  G4double pcm[3][18]; // pcm is an array with 3 rows and vecLen columns
330 
331  G4double totalMass = 0.0;
332  G4double extraMass = 0;
333  G4double sm[18];
334 
335  for (i=0; i<vecLen; ++i) {
336  mass[i] = vec[i]->GetMass()/GeV;
337  if(vec[i]->GetSide() == -2) extraMass+=vec[i]->GetMass()/GeV;
338  vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
339  pcm[0][i] = 0.0; // x-momentum of i-th particle
340  pcm[1][i] = 0.0; // y-momentum of i-th particle
341  pcm[2][i] = 0.0; // z-momentum of i-th particle
342  energy[i] = mass[i]; // total energy of i-th particle
343  totalMass += mass[i];
344  sm[i] = totalMass;
345  }
346 
347  G4double totalE = totalEnergy/GeV;
348  if (totalMass > totalE) {
349  //G4cerr << "*** Error in G4RPGReaction::GenerateNBodyEvent" << G4endl;
350  //G4cerr << " total mass (" << totalMass*GeV << "MeV) > total energy ("
351  // << totalEnergy << "MeV)" << G4endl;
352  totalE = totalMass;
353  return -1.0;
354  }
355 
356  G4double kineticEnergy = totalE - totalMass;
357  G4double emm[18];
358  emm[0] = mass[0];
359  emm[vecLen-1] = totalE;
360 
361  if (vecLen > 2) { // the random numbers are sorted
362  G4double ran[18];
363  for( i=0; i<vecLen; ++i )ran[i] = G4UniformRand();
364  for (i=0; i<vecLen-2; ++i) {
365  for (G4int j=vecLen-2; j>i; --j) {
366  if (ran[i] > ran[j]) {
367  G4double temp = ran[i];
368  ran[i] = ran[j];
369  ran[j] = temp;
370  }
371  }
372  }
373  for( i=1; i<vecLen-1; ++i )emm[i] = ran[i-1]*kineticEnergy + sm[i];
374  }
375 
376  // Weight is the sum of logarithms of terms instead of the product of terms
377 
378  G4bool lzero = true;
379  G4double wtmax = 0.0;
380  if (constantCrossSection) {
381  G4double emmax = kineticEnergy + mass[0];
382  G4double emmin = 0.0;
383  for( i=1; i<vecLen; ++i )
384  {
385  emmin += mass[i-1];
386  emmax += mass[i];
387  G4double wtfc = 0.0;
388  if( emmax*emmax > 0.0 )
389  {
390  G4double arg = emmax*emmax
391  + (emmin*emmin-mass[i]*mass[i])*(emmin*emmin-mass[i]*mass[i])/(emmax*emmax)
392  - 2.0*(emmin*emmin+mass[i]*mass[i]);
393  if( arg > 0.0 )wtfc = 0.5*std::sqrt( arg );
394  }
395  if( wtfc == 0.0 )
396  {
397  lzero = false;
398  break;
399  }
400  wtmax += std::log( wtfc );
401  }
402  if( lzero )
403  wtmax = -wtmax;
404  else
405  wtmax = expxu;
406  } else {
407  // ffq(n) = pi*(2*pi)^(n-2)/(n-2)!
408  const G4double ffq[18] = { 0., 3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
409  256.3704, 268.4705, 240.9780, 189.2637,
410  132.1308, 83.0202, 47.4210, 24.8295,
411  12.0006, 5.3858, 2.2560, 0.8859 };
412  wtmax = std::log( std::pow( kineticEnergy, vecLen-2 ) * ffq[vecLen-1] / totalE );
413  }
414 
415  // Calculate momenta for secondaries
416 
417  lzero = true;
418  G4double pd[50];
419 
420  for( i=0; i<vecLen-1; ++i )
421  {
422  pd[i] = 0.0;
423  if( emm[i+1]*emm[i+1] > 0.0 )
424  {
425  G4double arg = emm[i+1]*emm[i+1]
426  + (emm[i]*emm[i]-mass[i+1]*mass[i+1])*(emm[i]*emm[i]-mass[i+1]*mass[i+1])
427  /(emm[i+1]*emm[i+1])
428  - 2.0*(emm[i]*emm[i]+mass[i+1]*mass[i+1]);
429  if( arg > 0.0 )pd[i] = 0.5*std::sqrt( arg );
430  }
431  if( pd[i] <= 0.0 ) // changed from == on 02 April 98
432  lzero = false;
433  else
434  wtmax += std::log( pd[i] );
435  }
436  G4double weight = 0.0; // weight is returned by GenerateNBodyEvent
437  if( lzero )weight = std::exp( std::max(std::min(wtmax,expxu),expxl) );
438 
439  G4double bang, cb, sb, s0, s1, s2, c, esys, a, b, gama, beta;
440  G4double ss;
441  pcm[0][0] = 0.0;
442  pcm[1][0] = pd[0];
443  pcm[2][0] = 0.0;
444  for( i=1; i<vecLen; ++i )
445  {
446  pcm[0][i] = 0.0;
447  pcm[1][i] = -pd[i-1];
448  pcm[2][i] = 0.0;
449  bang = twopi*G4UniformRand();
450  cb = std::cos(bang);
451  sb = std::sin(bang);
452  c = 2.0*G4UniformRand() - 1.0;
453  ss = std::sqrt( std::fabs( 1.0-c*c ) );
454  if( i < vecLen-1 )
455  {
456  esys = std::sqrt(pd[i]*pd[i] + emm[i]*emm[i]);
457  beta = pd[i]/esys;
458  gama = esys/emm[i];
459  for( G4int j=0; j<=i; ++j )
460  {
461  s0 = pcm[0][j];
462  s1 = pcm[1][j];
463  s2 = pcm[2][j];
464  energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
465  a = s0*c - s1*ss; // rotation
466  pcm[1][j] = s0*ss + s1*c;
467  b = pcm[2][j];
468  pcm[0][j] = a*cb - b*sb;
469  pcm[2][j] = a*sb + b*cb;
470  pcm[1][j] = gama*(pcm[1][j] + beta*energy[j]);
471  }
472  }
473  else
474  {
475  for( G4int j=0; j<=i; ++j )
476  {
477  s0 = pcm[0][j];
478  s1 = pcm[1][j];
479  s2 = pcm[2][j];
480  energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
481  a = s0*c - s1*s; // rotation
482  pcm[1][j] = s0*ss + s1*c;
483  b = pcm[2][j];
484  pcm[0][j] = a*cb - b*sb;
485  pcm[2][j] = a*sb + b*cb;
486  }
487  }
488  }
489 
490  for (i=0; i<vecLen; ++i) {
491  vec[i]->SetMomentum( pcm[0][i]*GeV, pcm[1][i]*GeV, pcm[2][i]*GeV );
492  vec[i]->SetTotalEnergy( energy[i]*GeV );
493  }
494 
495  return weight;
496 }
497 
498 
499 G4double
500 G4RPGReaction::GenerateNBodyEventT(const G4double totalEnergy,
501  const G4bool constantCrossSection,
502  std::vector<G4ReactionProduct*>& tempList)
503 {
504  G4int i;
505  const G4double expxu = 82.; // upper bound for arg. of exp
506  const G4double expxl = -expxu; // lower bound for arg. of exp
507  G4int listLen = tempList.size();
508 
509  if (listLen < 2) {
510  G4cerr << "*** Error in G4RPGReaction::GenerateNBodyEvent" << G4endl;
511  G4cerr << " number of particles < 2" << G4endl;
512  G4cerr << "totalEnergy = " << totalEnergy << "MeV, listLen = " << listLen << G4endl;
513  return -1.0;
514  }
515 
516  G4double mass[18]; // mass of each particle
517  G4double energy[18]; // total energy of each particle
518  G4double pcm[3][18]; // pcm is an array with 3 rows and listLen columns
519 
520  G4double totalMass = 0.0;
521  G4double extraMass = 0;
522  G4double sm[18];
523 
524  for (i=0; i<listLen; ++i) {
525  mass[i] = tempList[i]->GetMass()/GeV;
526  if(tempList[i]->GetSide() == -2) extraMass+=tempList[i]->GetMass()/GeV;
527  tempList[i]->SetMomentum( 0.0, 0.0, 0.0 );
528  pcm[0][i] = 0.0; // x-momentum of i-th particle
529  pcm[1][i] = 0.0; // y-momentum of i-th particle
530  pcm[2][i] = 0.0; // z-momentum of i-th particle
531  energy[i] = mass[i]; // total energy of i-th particle
532  totalMass += mass[i];
533  sm[i] = totalMass;
534  }
535 
536  G4double totalE = totalEnergy/GeV;
537  if (totalMass > totalE) {
538  totalE = totalMass;
539  return -1.0;
540  }
541 
542  G4double kineticEnergy = totalE - totalMass;
543  G4double emm[18];
544  emm[0] = mass[0];
545  emm[listLen-1] = totalE;
546 
547  if (listLen > 2) { // the random numbers are sorted
548  G4double ran[18];
549  for( i=0; i<listLen; ++i )ran[i] = G4UniformRand();
550  for (i=0; i<listLen-2; ++i) {
551  for (G4int j=listLen-2; j>i; --j) {
552  if (ran[i] > ran[j]) {
553  G4double temp = ran[i];
554  ran[i] = ran[j];
555  ran[j] = temp;
556  }
557  }
558  }
559  for( i=1; i<listLen-1; ++i )emm[i] = ran[i-1]*kineticEnergy + sm[i];
560  }
561 
562  // Weight is the sum of logarithms of terms instead of the product of terms
563 
564  G4bool lzero = true;
565  G4double wtmax = 0.0;
566  if (constantCrossSection) {
567  G4double emmax = kineticEnergy + mass[0];
568  G4double emmin = 0.0;
569  for( i=1; i<listLen; ++i )
570  {
571  emmin += mass[i-1];
572  emmax += mass[i];
573  G4double wtfc = 0.0;
574  if( emmax*emmax > 0.0 )
575  {
576  G4double arg = emmax*emmax
577  + (emmin*emmin-mass[i]*mass[i])*(emmin*emmin-mass[i]*mass[i])/(emmax*emmax)
578  - 2.0*(emmin*emmin+mass[i]*mass[i]);
579  if( arg > 0.0 )wtfc = 0.5*std::sqrt( arg );
580  }
581  if( wtfc == 0.0 )
582  {
583  lzero = false;
584  break;
585  }
586  wtmax += std::log( wtfc );
587  }
588  if( lzero )
589  wtmax = -wtmax;
590  else
591  wtmax = expxu;
592  } else {
593  // ffq(n) = pi*(2*pi)^(n-2)/(n-2)!
594  const G4double ffq[18] = { 0., 3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
595  256.3704, 268.4705, 240.9780, 189.2637,
596  132.1308, 83.0202, 47.4210, 24.8295,
597  12.0006, 5.3858, 2.2560, 0.8859 };
598  wtmax = std::log( std::pow( kineticEnergy, listLen-2 ) * ffq[listLen-1] / totalE );
599  }
600 
601  // Calculate momenta for secondaries
602 
603  lzero = true;
604  G4double pd[50];
605 
606  for( i=0; i<listLen-1; ++i )
607  {
608  pd[i] = 0.0;
609  if( emm[i+1]*emm[i+1] > 0.0 )
610  {
611  G4double arg = emm[i+1]*emm[i+1]
612  + (emm[i]*emm[i]-mass[i+1]*mass[i+1])*(emm[i]*emm[i]-mass[i+1]*mass[i+1])
613  /(emm[i+1]*emm[i+1])
614  - 2.0*(emm[i]*emm[i]+mass[i+1]*mass[i+1]);
615  if( arg > 0.0 )pd[i] = 0.5*std::sqrt( arg );
616  }
617  if( pd[i] <= 0.0 ) // changed from == on 02 April 98
618  lzero = false;
619  else
620  wtmax += std::log( pd[i] );
621  }
622  G4double weight = 0.0; // weight is returned by GenerateNBodyEvent
623  if( lzero )weight = std::exp( std::max(std::min(wtmax,expxu),expxl) );
624 
625  G4double bang, cb, sb, s0, s1, s2, c, esys, a, b, gama, beta;
626  G4double ss;
627  pcm[0][0] = 0.0;
628  pcm[1][0] = pd[0];
629  pcm[2][0] = 0.0;
630  for( i=1; i<listLen; ++i )
631  {
632  pcm[0][i] = 0.0;
633  pcm[1][i] = -pd[i-1];
634  pcm[2][i] = 0.0;
635  bang = twopi*G4UniformRand();
636  cb = std::cos(bang);
637  sb = std::sin(bang);
638  c = 2.0*G4UniformRand() - 1.0;
639  ss = std::sqrt( std::fabs( 1.0-c*c ) );
640  if( i < listLen-1 )
641  {
642  esys = std::sqrt(pd[i]*pd[i] + emm[i]*emm[i]);
643  beta = pd[i]/esys;
644  gama = esys/emm[i];
645  for( G4int j=0; j<=i; ++j )
646  {
647  s0 = pcm[0][j];
648  s1 = pcm[1][j];
649  s2 = pcm[2][j];
650  energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
651  a = s0*c - s1*ss; // rotation
652  pcm[1][j] = s0*ss + s1*c;
653  b = pcm[2][j];
654  pcm[0][j] = a*cb - b*sb;
655  pcm[2][j] = a*sb + b*cb;
656  pcm[1][j] = gama*(pcm[1][j] + beta*energy[j]);
657  }
658  }
659  else
660  {
661  for( G4int j=0; j<=i; ++j )
662  {
663  s0 = pcm[0][j];
664  s1 = pcm[1][j];
665  s2 = pcm[2][j];
666  energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
667  a = s0*c - s1*ss; // rotation
668  pcm[1][j] = s0*ss + s1*c;
669  b = pcm[2][j];
670  pcm[0][j] = a*cb - b*sb;
671  pcm[2][j] = a*sb + b*cb;
672  }
673  }
674  }
675 
676  for (i=0; i<listLen; ++i) {
677  tempList[i]->SetMomentum(pcm[0][i]*GeV, pcm[1][i]*GeV, pcm[2][i]*GeV);
678  tempList[i]->SetTotalEnergy(energy[i]*GeV);
679  }
680 
681  return weight;
682 }
683 
684 
685 G4double G4RPGReaction::normal()
686 {
687  G4double ran = -6.0;
688  for( G4int i=0; i<12; ++i )ran += G4UniformRand();
689  return ran;
690 }
691 
692 
693 void G4RPGReaction::Defs1(const G4ReactionProduct& modifiedOriginal,
694  G4ReactionProduct& currentParticle,
695  G4ReactionProduct& targetParticle,
696  G4FastVector<G4ReactionProduct,256>& vec,
697  G4int& vecLen)
698 {
699  // Rotate final state particle momenta by initial particle direction
700 
701  const G4double pjx = modifiedOriginal.GetMomentum().x()/MeV;
702  const G4double pjy = modifiedOriginal.GetMomentum().y()/MeV;
703  const G4double pjz = modifiedOriginal.GetMomentum().z()/MeV;
704  const G4double p = modifiedOriginal.GetMomentum().mag()/MeV;
705 
706  if (pjx*pjx+pjy*pjy > 0.0) {
707  G4double cost, sint, ph, cosp, sinp, pix, piy, piz;
708  cost = pjz/p;
709  sint = std::sqrt(std::abs((1.0-cost)*(1.0+cost)));
710  if( pjy < 0.0 )
711  ph = 3*halfpi;
712  else
713  ph = halfpi;
714  if( std::abs( pjx ) > 0.001*MeV )ph = std::atan2(pjy,pjx);
715  cosp = std::cos(ph);
716  sinp = std::sin(ph);
717  pix = currentParticle.GetMomentum().x()/MeV;
718  piy = currentParticle.GetMomentum().y()/MeV;
719  piz = currentParticle.GetMomentum().z()/MeV;
720  currentParticle.SetMomentum((cost*cosp*pix - sinp*piy + sint*cosp*piz)*MeV,
721  (cost*sinp*pix + cosp*piy + sint*sinp*piz)*MeV,
722  (-sint*pix + cost*piz)*MeV);
723  pix = targetParticle.GetMomentum().x()/MeV;
724  piy = targetParticle.GetMomentum().y()/MeV;
725  piz = targetParticle.GetMomentum().z()/MeV;
726  targetParticle.SetMomentum((cost*cosp*pix - sinp*piy + sint*cosp*piz)*MeV,
727  (cost*sinp*pix + cosp*piy + sint*sinp*piz)*MeV,
728  (-sint*pix + cost*piz)*MeV);
729 
730  for (G4int i=0; i<vecLen; ++i) {
731  pix = vec[i]->GetMomentum().x()/MeV;
732  piy = vec[i]->GetMomentum().y()/MeV;
733  piz = vec[i]->GetMomentum().z()/MeV;
734  vec[i]->SetMomentum((cost*cosp*pix - sinp*piy + sint*cosp*piz)*MeV,
735  (cost*sinp*pix + cosp*piy + sint*sinp*piz)*MeV,
736  (-sint*pix + cost*piz)*MeV);
737  }
738 
739  } else {
740  if (pjz < 0.0) {
741  currentParticle.SetMomentum( -currentParticle.GetMomentum().z() );
742  targetParticle.SetMomentum( -targetParticle.GetMomentum().z() );
743  for (G4int i=0; i<vecLen; ++i) vec[i]->SetMomentum( -vec[i]->GetMomentum().z() );
744  }
745  }
746 }
747 
748 
749  void G4RPGReaction::Rotate(
750  const G4double numberofFinalStateNucleons,
751  const G4ThreeVector &temp,
752  const G4ReactionProduct &modifiedOriginal, // Fermi motion & evap. effect included
753  const G4HadProjectile *originalIncident, // original incident particle
754  const G4Nucleus &targetNucleus,
755  G4ReactionProduct &currentParticle,
756  G4ReactionProduct &targetParticle,
757  G4FastVector<G4ReactionProduct,256> &vec,
758  G4int &vecLen )
759  {
760  // derived from original FORTRAN code in GENXPT and TWOCLU by H. Fesefeldt
761  //
762  // Rotate in direction of z-axis, this does disturb in some way our
763  // inclusive distributions, but it is necessary for momentum conservation
764  //
765  const G4double atomicWeight = targetNucleus.GetA_asInt();
766  const G4double logWeight = std::log(atomicWeight);
767 
768  G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
769  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
770  G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
771 
772  G4int i;
773  G4ThreeVector pseudoParticle[4];
774  for( i=0; i<4; ++i )pseudoParticle[i].set(0,0,0);
775  pseudoParticle[0] = currentParticle.GetMomentum()
776  + targetParticle.GetMomentum();
777  for( i=0; i<vecLen; ++i )
778  pseudoParticle[0] = pseudoParticle[0] + (vec[i]->GetMomentum());
779  //
780  // Some smearing in transverse direction from Fermi motion
781  //
782  G4float pp, pp1;
783  G4double alekw, p;
784  G4double r1, r2, a1, ran1, ran2, xxh, exh, pxTemp, pyTemp, pzTemp;
785 
786  r1 = twopi*G4UniformRand();
787  r2 = G4UniformRand();
788  a1 = std::sqrt(-2.0*std::log(r2));
789  ran1 = a1*std::sin(r1)*0.020*numberofFinalStateNucleons*GeV;
790  ran2 = a1*std::cos(r1)*0.020*numberofFinalStateNucleons*GeV;
791  G4ThreeVector fermir(ran1, ran2, 0);
792 
793  pseudoParticle[0] = pseudoParticle[0]+fermir; // all particles + fermir
794  pseudoParticle[2] = temp; // original in cms system
795  pseudoParticle[3] = pseudoParticle[0];
796 
797  pseudoParticle[1] = pseudoParticle[2].cross(pseudoParticle[3]);
798  G4double rotation = 2.*pi*G4UniformRand();
799  pseudoParticle[1] = pseudoParticle[1].rotate(rotation, pseudoParticle[3]);
800  pseudoParticle[2] = pseudoParticle[3].cross(pseudoParticle[1]);
801  for(G4int ii=1; ii<=3; ii++)
802  {
803  p = pseudoParticle[ii].mag();
804  if( p == 0.0 )
805  pseudoParticle[ii]= G4ThreeVector( 0.0, 0.0, 0.0 );
806  else
807  pseudoParticle[ii]= pseudoParticle[ii] * (1./p);
808  }
809 
810  pxTemp = pseudoParticle[1].dot(currentParticle.GetMomentum());
811  pyTemp = pseudoParticle[2].dot(currentParticle.GetMomentum());
812  pzTemp = pseudoParticle[3].dot(currentParticle.GetMomentum());
813  currentParticle.SetMomentum( pxTemp, pyTemp, pzTemp );
814 
815  pxTemp = pseudoParticle[1].dot(targetParticle.GetMomentum());
816  pyTemp = pseudoParticle[2].dot(targetParticle.GetMomentum());
817  pzTemp = pseudoParticle[3].dot(targetParticle.GetMomentum());
818  targetParticle.SetMomentum( pxTemp, pyTemp, pzTemp );
819 
820  for( i=0; i<vecLen; ++i )
821  {
822  pxTemp = pseudoParticle[1].dot(vec[i]->GetMomentum());
823  pyTemp = pseudoParticle[2].dot(vec[i]->GetMomentum());
824  pzTemp = pseudoParticle[3].dot(vec[i]->GetMomentum());
825  vec[i]->SetMomentum( pxTemp, pyTemp, pzTemp );
826  }
827  //
828  // Rotate in direction of primary particle, subtract binding energies
829  // and make some further corrections if required
830  //
831  Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
832  G4double ekin;
833  G4double dekin = 0.0;
834  G4double ek1 = 0.0;
835  G4int npions = 0;
836  if( atomicWeight >= 1.5 ) // self-absorption in heavy molecules
837  {
838  // corrections for single particle spectra (shower particles)
839  //
840  const G4double alem[] = { 1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00 };
841  const G4double val0[] = { 0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65 };
842  alekw = std::log( originalIncident->GetKineticEnergy()/GeV );
843  exh = 1.0;
844  if( alekw > alem[0] ) // get energy bin
845  {
846  exh = val0[6];
847  for( G4int j=1; j<7; ++j )
848  {
849  if( alekw < alem[j] ) // use linear interpolation/extrapolation
850  {
851  G4double rcnve = (val0[j] - val0[j-1]) / (alem[j] - alem[j-1]);
852  exh = rcnve * alekw + val0[j-1] - rcnve * alem[j-1];
853  break;
854  }
855  }
856  exh = 1.0 - exh;
857  }
858  const G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
859  ekin = currentParticle.GetKineticEnergy()/GeV - cfa*(1+normal()/2.0);
860  ekin = std::max( 1.0e-6, ekin );
861  xxh = 1.0;
862  if( (modifiedOriginal.GetDefinition() == aPiPlus ||
863  modifiedOriginal.GetDefinition() == aPiMinus) &&
864  currentParticle.GetDefinition() == aPiZero &&
865  G4UniformRand() <= logWeight) xxh = exh;
866  dekin += ekin*(1.0-xxh);
867  ekin *= xxh;
868  if (currentParticle.GetDefinition()->GetParticleSubType() == "pi") {
869  ++npions;
870  ek1 += ekin;
871  }
872  currentParticle.SetKineticEnergy( ekin*GeV );
873  pp = currentParticle.GetTotalMomentum()/MeV;
874  pp1 = currentParticle.GetMomentum().mag()/MeV;
875  if( pp1 < 0.001*MeV )
876  {
877  G4double costheta = 2.*G4UniformRand() - 1.;
878  G4double sintheta = std::sqrt(1. - costheta*costheta);
879  G4double phi = twopi*G4UniformRand();
880  currentParticle.SetMomentum( pp*sintheta*std::cos(phi)*MeV,
881  pp*sintheta*std::sin(phi)*MeV,
882  pp*costheta*MeV ) ;
883  }
884  else
885  currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
886  ekin = targetParticle.GetKineticEnergy()/GeV - cfa*(1+normal()/2.0);
887  ekin = std::max( 1.0e-6, ekin );
888  xxh = 1.0;
889  if( (modifiedOriginal.GetDefinition() == aPiPlus ||
890  modifiedOriginal.GetDefinition() == aPiMinus) &&
891  targetParticle.GetDefinition() == aPiZero &&
892  G4UniformRand() < logWeight) xxh = exh;
893  dekin += ekin*(1.0-xxh);
894  ekin *= xxh;
895  if (targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
896  ++npions;
897  ek1 += ekin;
898  }
899  targetParticle.SetKineticEnergy( ekin*GeV );
900  pp = targetParticle.GetTotalMomentum()/MeV;
901  pp1 = targetParticle.GetMomentum().mag()/MeV;
902  if( pp1 < 0.001*MeV )
903  {
904  G4double costheta = 2.*G4UniformRand() - 1.;
905  G4double sintheta = std::sqrt(1. - costheta*costheta);
906  G4double phi = twopi*G4UniformRand();
907  targetParticle.SetMomentum( pp*sintheta*std::cos(phi)*MeV,
908  pp*sintheta*std::sin(phi)*MeV,
909  pp*costheta*MeV ) ;
910  }
911  else
912  targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
913  for( i=0; i<vecLen; ++i )
914  {
915  ekin = vec[i]->GetKineticEnergy()/GeV - cfa*(1+normal()/2.0);
916  ekin = std::max( 1.0e-6, ekin );
917  xxh = 1.0;
918  if( (modifiedOriginal.GetDefinition() == aPiPlus ||
919  modifiedOriginal.GetDefinition() == aPiMinus) &&
920  vec[i]->GetDefinition() == aPiZero &&
921  G4UniformRand() < logWeight) xxh = exh;
922  dekin += ekin*(1.0-xxh);
923  ekin *= xxh;
924  if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
925  ++npions;
926  ek1 += ekin;
927  }
928  vec[i]->SetKineticEnergy( ekin*GeV );
929  pp = vec[i]->GetTotalMomentum()/MeV;
930  pp1 = vec[i]->GetMomentum().mag()/MeV;
931  if( pp1 < 0.001*MeV )
932  {
933  G4double costheta = 2.*G4UniformRand() - 1.;
934  G4double sintheta = std::sqrt(1. - costheta*costheta);
935  G4double phi = twopi*G4UniformRand();
936  vec[i]->SetMomentum( pp*sintheta*std::cos(phi)*MeV,
937  pp*sintheta*std::sin(phi)*MeV,
938  pp*costheta*MeV ) ;
939  }
940  else
941  vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
942  }
943  }
944  if( (ek1 != 0.0) && (npions > 0) )
945  {
946  dekin = 1.0 + dekin/ek1;
947  //
948  // first do the incident particle
949  //
950  if (currentParticle.GetDefinition()->GetParticleSubType() == "pi")
951  {
952  currentParticle.SetKineticEnergy(
953  std::max( 0.001*MeV, dekin*currentParticle.GetKineticEnergy() ) );
954  pp = currentParticle.GetTotalMomentum()/MeV;
955  pp1 = currentParticle.GetMomentum().mag()/MeV;
956  if( pp1 < 0.001 )
957  {
958  G4double costheta = 2.*G4UniformRand() - 1.;
959  G4double sintheta = std::sqrt(1. - costheta*costheta);
960  G4double phi = twopi*G4UniformRand();
961  currentParticle.SetMomentum( pp*sintheta*std::cos(phi)*MeV,
962  pp*sintheta*std::sin(phi)*MeV,
963  pp*costheta*MeV ) ;
964  } else {
965  currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
966  }
967  }
968 
969  if (targetParticle.GetDefinition()->GetParticleSubType() == "pi")
970  {
971  targetParticle.SetKineticEnergy(
972  std::max( 0.001*MeV, dekin*targetParticle.GetKineticEnergy() ) );
973  pp = targetParticle.GetTotalMomentum()/MeV;
974  pp1 = targetParticle.GetMomentum().mag()/MeV;
975  if( pp1 < 0.001 )
976  {
977  G4double costheta = 2.*G4UniformRand() - 1.;
978  G4double sintheta = std::sqrt(1. - costheta*costheta);
979  G4double phi = twopi*G4UniformRand();
980  targetParticle.SetMomentum( pp*sintheta*std::cos(phi)*MeV,
981  pp*sintheta*std::sin(phi)*MeV,
982  pp*costheta*MeV ) ;
983  } else {
984  targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
985  }
986  }
987 
988  for( i=0; i<vecLen; ++i )
989  {
990  if (vec[i]->GetDefinition()->GetParticleSubType() == "pi")
991  {
992  vec[i]->SetKineticEnergy( std::max( 0.001*MeV, dekin*vec[i]->GetKineticEnergy() ) );
993  pp = vec[i]->GetTotalMomentum()/MeV;
994  pp1 = vec[i]->GetMomentum().mag()/MeV;
995  if( pp1 < 0.001 )
996  {
997  G4double costheta = 2.*G4UniformRand() - 1.;
998  G4double sintheta = std::sqrt(1. - costheta*costheta);
999  G4double phi = twopi*G4UniformRand();
1000  vec[i]->SetMomentum( pp*sintheta*std::cos(phi)*MeV,
1001  pp*sintheta*std::sin(phi)*MeV,
1002  pp*costheta*MeV ) ;
1003  } else {
1004  vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
1005  }
1006  }
1007  } // for i
1008  } // if (ek1 != 0)
1009  }
1010 
1011  std::pair<G4int, G4int> G4RPGReaction::GetFinalStateNucleons(
1012  const G4DynamicParticle* originalTarget,
1013  const G4FastVector<G4ReactionProduct,256>& vec,
1014  const G4int& vecLen)
1015  {
1016  // Get number of protons and neutrons removed from the target nucleus
1017 
1018  G4int protonsRemoved = 0;
1019  G4int neutronsRemoved = 0;
1020  if (originalTarget->GetDefinition()->GetParticleName() == "proton")
1021  protonsRemoved++;
1022  else
1023  neutronsRemoved++;
1024 
1025  G4String secName;
1026  for (G4int i = 0; i < vecLen; i++) {
1027  secName = vec[i]->GetDefinition()->GetParticleName();
1028  if (secName == "proton") {
1029  protonsRemoved++;
1030  } else if (secName == "neutron") {
1031  neutronsRemoved++;
1032  } else if (secName == "anti_proton") {
1033  protonsRemoved--;
1034  } else if (secName == "anti_neutron") {
1035  neutronsRemoved--;
1036  }
1037  }
1038 
1039  return std::pair<G4int, G4int>(protonsRemoved, neutronsRemoved);
1040  }
1041 
1042 
1043  G4ThreeVector G4RPGReaction::Isotropic(const G4double& pp)
1044  {
1045  G4double costheta = 2.*G4UniformRand() - 1.;
1046  G4double sintheta = std::sqrt(1. - costheta*costheta);
1047  G4double phi = twopi*G4UniformRand();
1048  return G4ThreeVector(pp*sintheta*std::cos(phi),
1049  pp*sintheta*std::sin(phi),
1050  pp*costheta);
1051  }
1052 
1053 
1054  void G4RPGReaction::MomentumCheck(
1055  const G4ReactionProduct &modifiedOriginal,
1056  G4ReactionProduct &currentParticle,
1057  G4ReactionProduct &targetParticle,
1058  G4FastVector<G4ReactionProduct,256> &vec,
1059  G4int &vecLen )
1060  {
1061  const G4double pOriginal = modifiedOriginal.GetTotalMomentum()/MeV;
1062  G4double testMomentum = currentParticle.GetMomentum().mag()/MeV;
1063  G4double pMass;
1064  if( testMomentum >= pOriginal )
1065  {
1066  pMass = currentParticle.GetMass()/MeV;
1067  currentParticle.SetTotalEnergy(
1068  std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
1069  currentParticle.SetMomentum(
1070  currentParticle.GetMomentum() * (pOriginal/testMomentum) );
1071  }
1072  testMomentum = targetParticle.GetMomentum().mag()/MeV;
1073  if( testMomentum >= pOriginal )
1074  {
1075  pMass = targetParticle.GetMass()/MeV;
1076  targetParticle.SetTotalEnergy(
1077  std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
1078  targetParticle.SetMomentum(
1079  targetParticle.GetMomentum() * (pOriginal/testMomentum) );
1080  }
1081  for( G4int i=0; i<vecLen; ++i )
1082  {
1083  testMomentum = vec[i]->GetMomentum().mag()/MeV;
1084  if( testMomentum >= pOriginal )
1085  {
1086  pMass = vec[i]->GetMass()/MeV;
1087  vec[i]->SetTotalEnergy(
1088  std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
1089  vec[i]->SetMomentum( vec[i]->GetMomentum() * (pOriginal/testMomentum) );
1090  }
1091  }
1092  }
1093 
1094  void G4RPGReaction::NuclearReaction(
1095  G4FastVector<G4ReactionProduct,4> &vec,
1096  G4int &vecLen,
1097  const G4HadProjectile *originalIncident,
1098  const G4Nucleus &targetNucleus,
1099  const G4double theAtomicMass,
1100  const G4double *mass )
1101  {
1102  // derived from original FORTRAN code NUCREC by H. Fesefeldt (12-Feb-1987)
1103  //
1104  // Nuclear reaction kinematics at low energies
1105  //
1106  G4ParticleDefinition *aGamma = G4Gamma::Gamma();
1107  G4ParticleDefinition *aProton = G4Proton::Proton();
1108  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
1109  G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
1110  G4ParticleDefinition *aTriton = G4Triton::Triton();
1111  G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
1112 
1113  const G4double aProtonMass = aProton->GetPDGMass()/MeV;
1114  const G4double aNeutronMass = aNeutron->GetPDGMass()/MeV;
1115  const G4double aDeuteronMass = aDeuteron->GetPDGMass()/MeV;
1116  const G4double aTritonMass = aTriton->GetPDGMass()/MeV;
1117  const G4double anAlphaMass = anAlpha->GetPDGMass()/MeV;
1118 
1119  G4ReactionProduct currentParticle;
1120  currentParticle = *originalIncident;
1121  //
1122  // Set beam particle, take kinetic energy of current particle as the
1123  // fundamental quantity. Due to the difficult kinematic, all masses have to
1124  // be assigned the best measured values
1125  //
1126  G4double p = currentParticle.GetTotalMomentum();
1127  G4double pp = currentParticle.GetMomentum().mag();
1128  if( pp <= 0.001*MeV )
1129  {
1130  G4double phinve = twopi*G4UniformRand();
1131  G4double rthnve = std::acos( std::max( -1.0, std::min( 1.0, -1.0 + 2.0*G4UniformRand() ) ) );
1132  currentParticle.SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
1133  p*std::sin(rthnve)*std::sin(phinve),
1134  p*std::cos(rthnve) );
1135  }
1136  else
1137  currentParticle.SetMomentum( currentParticle.GetMomentum() * (p/pp) );
1138  //
1139  // calculate Q-value of reactions
1140  //
1141  G4double currentKinetic = currentParticle.GetKineticEnergy()/MeV;
1142  G4double currentMass = currentParticle.GetDefinition()->GetPDGMass()/MeV;
1143  G4double qv = currentKinetic + theAtomicMass + currentMass;
1144 
1145  G4double qval[9];
1146  qval[0] = qv - mass[0];
1147  qval[1] = qv - mass[1] - aNeutronMass;
1148  qval[2] = qv - mass[2] - aProtonMass;
1149  qval[3] = qv - mass[3] - aDeuteronMass;
1150  qval[4] = qv - mass[4] - aTritonMass;
1151  qval[5] = qv - mass[5] - anAlphaMass;
1152  qval[6] = qv - mass[6] - aNeutronMass - aNeutronMass;
1153  qval[7] = qv - mass[7] - aNeutronMass - aProtonMass;
1154  qval[8] = qv - mass[8] - aProtonMass - aProtonMass;
1155 
1156  if( currentParticle.GetDefinition() == aNeutron )
1157  {
1158  const G4double A = targetNucleus.GetA_asInt(); // atomic weight
1159  if( G4UniformRand() > ((A-1.0)/230.0)*((A-1.0)/230.0) )
1160  qval[0] = 0.0;
1161  if( G4UniformRand() >= currentKinetic/7.9254*A )
1162  qval[2] = qval[3] = qval[4] = qval[5] = qval[8] = 0.0;
1163  }
1164  else
1165  qval[0] = 0.0;
1166 
1167  G4int i;
1168  qv = 0.0;
1169  for( i=0; i<9; ++i )
1170  {
1171  if( mass[i] < 500.0*MeV )qval[i] = 0.0;
1172  if( qval[i] < 0.0 )qval[i] = 0.0;
1173  qv += qval[i];
1174  }
1175  G4double qv1 = 0.0;
1176  G4double ran = G4UniformRand();
1177  G4int index;
1178  for( index=0; index<9; ++index )
1179  {
1180  if( qval[index] > 0.0 )
1181  {
1182  qv1 += qval[index]/qv;
1183  if( ran <= qv1 )break;
1184  }
1185  }
1186  if( index == 9 ) // loop continued to the end
1187  {
1188  throw G4HadronicException(__FILE__, __LINE__,
1189  "G4RPGReaction::NuclearReaction: inelastic reaction kinematically not possible");
1190  }
1191  G4double ke = currentParticle.GetKineticEnergy()/GeV;
1192  G4int nt = 2;
1193  if( (index>=6) || (G4UniformRand()<std::min(0.5,ke*10.0)) )nt = 3;
1194 
1195  G4ReactionProduct **v = new G4ReactionProduct * [3];
1196  v[0] = new G4ReactionProduct;
1197  v[1] = new G4ReactionProduct;
1198  v[2] = new G4ReactionProduct;
1199 
1200  v[0]->SetMass( mass[index]*MeV );
1201  switch( index )
1202  {
1203  case 0:
1204  v[1]->SetDefinition( aGamma );
1205  v[2]->SetDefinition( aGamma );
1206  break;
1207  case 1:
1208  v[1]->SetDefinition( aNeutron );
1209  v[2]->SetDefinition( aGamma );
1210  break;
1211  case 2:
1212  v[1]->SetDefinition( aProton );
1213  v[2]->SetDefinition( aGamma );
1214  break;
1215  case 3:
1216  v[1]->SetDefinition( aDeuteron );
1217  v[2]->SetDefinition( aGamma );
1218  break;
1219  case 4:
1220  v[1]->SetDefinition( aTriton );
1221  v[2]->SetDefinition( aGamma );
1222  break;
1223  case 5:
1224  v[1]->SetDefinition( anAlpha );
1225  v[2]->SetDefinition( aGamma );
1226  break;
1227  case 6:
1228  v[1]->SetDefinition( aNeutron );
1229  v[2]->SetDefinition( aNeutron );
1230  break;
1231  case 7:
1232  v[1]->SetDefinition( aNeutron );
1233  v[2]->SetDefinition( aProton );
1234  break;
1235  case 8:
1236  v[1]->SetDefinition( aProton );
1237  v[2]->SetDefinition( aProton );
1238  break;
1239  }
1240  //
1241  // calculate centre of mass energy
1242  //
1243  G4ReactionProduct pseudo1;
1244  pseudo1.SetMass( theAtomicMass*MeV );
1245  pseudo1.SetTotalEnergy( theAtomicMass*MeV );
1246  G4ReactionProduct pseudo2 = currentParticle + pseudo1;
1247  pseudo2.SetMomentum( pseudo2.GetMomentum() * (-1.0) );
1248  //
1249  // use phase space routine in centre of mass system
1250  //
1251  G4FastVector<G4ReactionProduct,256> tempV;
1252  tempV.Initialize( nt );
1253  G4int tempLen = 0;
1254  tempV.SetElement( tempLen++, v[0] );
1255  tempV.SetElement( tempLen++, v[1] );
1256  if( nt == 3 )tempV.SetElement( tempLen++, v[2] );
1257  G4bool constantCrossSection = true;
1258  GenerateNBodyEvent( pseudo2.GetMass()/MeV, constantCrossSection, tempV, tempLen );
1259  v[0]->Lorentz( *v[0], pseudo2 );
1260  v[1]->Lorentz( *v[1], pseudo2 );
1261  if( nt == 3 )v[2]->Lorentz( *v[2], pseudo2 );
1262 
1263  G4bool particleIsDefined = false;
1264  if( v[0]->GetMass()/MeV - aProtonMass < 0.1 )
1265  {
1266  v[0]->SetDefinition( aProton );
1267  particleIsDefined = true;
1268  }
1269  else if( v[0]->GetMass()/MeV - aNeutronMass < 0.1 )
1270  {
1271  v[0]->SetDefinition( aNeutron );
1272  particleIsDefined = true;
1273  }
1274  else if( v[0]->GetMass()/MeV - aDeuteronMass < 0.1 )
1275  {
1276  v[0]->SetDefinition( aDeuteron );
1277  particleIsDefined = true;
1278  }
1279  else if( v[0]->GetMass()/MeV - aTritonMass < 0.1 )
1280  {
1281  v[0]->SetDefinition( aTriton );
1282  particleIsDefined = true;
1283  }
1284  else if( v[0]->GetMass()/MeV - anAlphaMass < 0.1 )
1285  {
1286  v[0]->SetDefinition( anAlpha );
1287  particleIsDefined = true;
1288  }
1289  currentParticle.SetKineticEnergy(
1290  std::max( 0.001, currentParticle.GetKineticEnergy()/MeV ) );
1291  p = currentParticle.GetTotalMomentum();
1292  pp = currentParticle.GetMomentum().mag();
1293  if( pp <= 0.001*MeV )
1294  {
1295  G4double phinve = twopi*G4UniformRand();
1296  G4double rthnve = std::acos( std::max( -1.0, std::min( 1.0, -1.0 + 2.0*G4UniformRand() ) ) );
1297  currentParticle.SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
1298  p*std::sin(rthnve)*std::sin(phinve),
1299  p*std::cos(rthnve) );
1300  }
1301  else
1302  currentParticle.SetMomentum( currentParticle.GetMomentum() * (p/pp) );
1303 
1304  if( particleIsDefined )
1305  {
1306  v[0]->SetKineticEnergy(
1307  std::max( 0.001, 0.5*G4UniformRand()*v[0]->GetKineticEnergy()/MeV ) );
1308  p = v[0]->GetTotalMomentum();
1309  pp = v[0]->GetMomentum().mag();
1310  if( pp <= 0.001*MeV )
1311  {
1312  G4double phinve = twopi*G4UniformRand();
1313  G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
1314  v[0]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
1315  p*std::sin(rthnve)*std::sin(phinve),
1316  p*std::cos(rthnve) );
1317  }
1318  else
1319  v[0]->SetMomentum( v[0]->GetMomentum() * (p/pp) );
1320  }
1321  if( (v[1]->GetDefinition() == aDeuteron) ||
1322  (v[1]->GetDefinition() == aTriton) ||
1323  (v[1]->GetDefinition() == anAlpha) )
1324  v[1]->SetKineticEnergy(
1325  std::max( 0.001, 0.5*G4UniformRand()*v[1]->GetKineticEnergy()/MeV ) );
1326  else
1327  v[1]->SetKineticEnergy( std::max( 0.001, v[1]->GetKineticEnergy()/MeV ) );
1328 
1329  p = v[1]->GetTotalMomentum();
1330  pp = v[1]->GetMomentum().mag();
1331  if( pp <= 0.001*MeV )
1332  {
1333  G4double phinve = twopi*G4UniformRand();
1334  G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
1335  v[1]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
1336  p*std::sin(rthnve)*std::sin(phinve),
1337  p*std::cos(rthnve) );
1338  }
1339  else
1340  v[1]->SetMomentum( v[1]->GetMomentum() * (p/pp) );
1341 
1342  if( nt == 3 )
1343  {
1344  if( (v[2]->GetDefinition() == aDeuteron) ||
1345  (v[2]->GetDefinition() == aTriton) ||
1346  (v[2]->GetDefinition() == anAlpha) )
1347  v[2]->SetKineticEnergy(
1348  std::max( 0.001, 0.5*G4UniformRand()*v[2]->GetKineticEnergy()/MeV ) );
1349  else
1350  v[2]->SetKineticEnergy( std::max( 0.001, v[2]->GetKineticEnergy()/MeV ) );
1351 
1352  p = v[2]->GetTotalMomentum();
1353  pp = v[2]->GetMomentum().mag();
1354  if( pp <= 0.001*MeV )
1355  {
1356  G4double phinve = twopi*G4UniformRand();
1357  G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
1358  v[2]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
1359  p*std::sin(rthnve)*std::sin(phinve),
1360  p*std::cos(rthnve) );
1361  }
1362  else
1363  v[2]->SetMomentum( v[2]->GetMomentum() * (p/pp) );
1364  }
1365  G4int del;
1366  for(del=0; del<vecLen; del++) delete vec[del];
1367  vecLen = 0;
1368  if( particleIsDefined )
1369  {
1370  vec.SetElement( vecLen++, v[0] );
1371  }
1372  else
1373  {
1374  delete v[0];
1375  }
1376  vec.SetElement( vecLen++, v[1] );
1377  if( nt == 3 )
1378  {
1379  vec.SetElement( vecLen++, v[2] );
1380  }
1381  else
1382  {
1383  delete v[2];
1384  }
1385  delete [] v;
1386  return;
1387  }
1388 
1389  /* end of file */
1390 
Generated on Sat May 25 2013 14:34:08 for Geant4 by   doxygen 1.8.4