ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGReaction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGReaction.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 
28 #include <iostream>
29 
30 #include "G4RPGReaction.hh"
31 #include "G4Log.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*G4Exp((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 += G4Log( 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 = G4Log( 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 += G4Log( pd[i] );
435  }
436  G4double weight = 0.0; // weight is returned by GenerateNBodyEvent
437  if( lzero )weight = G4Exp( 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 += G4Log( 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 = G4Log( 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 += G4Log( pd[i] );
621  }
622  G4double weight = 0.0; // weight is returned by GenerateNBodyEvent
623  if( lzero )weight = G4Exp( 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 = G4Log(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  G4double 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*G4Log(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 = G4Log( 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.)*G4Exp(-(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