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