source: trunk/source/processes/hadronic/models/rpg/src/G4RPGReaction.cc@ 1199

Last change on this file since 1199 was 1196, checked in by garnier, 16 years ago

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