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

Last change on this file since 880 was 819, checked in by garnier, 17 years ago

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