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