Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | Protected Member Functions | List of all members
G4RPGReaction Class Reference

#include <G4RPGReaction.hh>

+ Inheritance diagram for G4RPGReaction:

Public Member Functions

 G4RPGReaction ()
 
virtual ~G4RPGReaction ()
 
G4bool ReactionStage (const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)
 
void AddBlackTrackParticles (const G4double, const G4int, const G4double, const G4int, const G4ReactionProduct &, G4int, G4int, const G4Nucleus &, G4FastVector< G4ReactionProduct, 256 > &, G4int &)
 
G4double GenerateNBodyEvent (const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
 
G4double GenerateNBodyEventT (const G4double totalEnergy, const G4bool constantCrossSection, std::vector< G4ReactionProduct * > &list)
 
void NuclearReaction (G4FastVector< G4ReactionProduct, 4 > &vec, G4int &vecLen, const G4HadProjectile *originalIncident, const G4Nucleus &aNucleus, const G4double theAtomicMass, const G4double *massVec)
 

Protected Member Functions

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)
 
void Defs1 (const G4ReactionProduct &modifiedOriginal, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
 
std::pair< G4int, G4int > GetFinalStateNucleons (const G4DynamicParticle *originalTarget, const G4FastVector< G4ReactionProduct, 256 > &vec, const G4int &vecLen)
 
void MomentumCheck (const G4ReactionProduct &modifiedOriginal, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
 
G4double normal ()
 
G4ThreeVector Isotropic (const G4double &)
 

Detailed Description

Definition at line 47 of file G4RPGReaction.hh.

Constructor & Destructor Documentation

◆ G4RPGReaction()

G4RPGReaction::G4RPGReaction ( )
inline

Definition at line 51 of file G4RPGReaction.hh.

51{}

◆ ~G4RPGReaction()

virtual G4RPGReaction::~G4RPGReaction ( )
inlinevirtual

Definition at line 53 of file G4RPGReaction.hh.

53{}

Member Function Documentation

◆ AddBlackTrackParticles()

void G4RPGReaction::AddBlackTrackParticles ( const G4double  epnb,
const G4int  npnb,
const G4double  edta,
const G4int  ndta,
const G4ReactionProduct &  modifiedOriginal,
G4int  PinNucleus,
G4int  NinNucleus,
const G4Nucleus &  targetNucleus,
G4FastVector< G4ReactionProduct, 256 > &  vec,
G4int &  vecLen 
)

Definition at line 56 of file G4RPGReaction.cc.

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}
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4FastVector::SetElement
void SetElement(G4int anIndex, Type *anElement)
Definition: G4FastVector.hh:76
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
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4ReactionProduct::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4ReactionProduct.hh:135
G4ReactionProduct::SetDefinition
void SetDefinition(G4ParticleDefinition *aParticleDefinition)
Definition: G4ReactionProduct.cc:155
G4Triton::Triton
static G4Triton * Triton()
Definition: G4Triton.cc:95

Referenced by G4RPGFragmentation::ReactionStage(), G4RPGTwoBody::ReactionStage(), and G4RPGTwoCluster::ReactionStage().

◆ Defs1()

void G4RPGReaction::Defs1 ( const G4ReactionProduct &  modifiedOriginal,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
G4FastVector< G4ReactionProduct, 256 > &  vec,
G4int &  vecLen 
)
protected

Definition at line 693 of file G4RPGReaction.cc.

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}
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::mag
double mag() const
G4ReactionProduct::SetMomentum
void SetMomentum(const G4double x, const G4double y, const G4double z)
Definition: G4ReactionProduct.cc:166
G4ReactionProduct::GetMomentum
G4ThreeVector GetMomentum() const
Definition: G4ReactionProduct.hh:120

Referenced by G4RPGTwoBody::ReactionStage(), and Rotate().

◆ GenerateNBodyEvent()

G4double G4RPGReaction::GenerateNBodyEvent ( const G4double  totalEnergy,
const G4bool  constantCrossSection,
G4FastVector< G4ReactionProduct, 256 > &  vec,
G4int &  vecLen 
)

Definition at line 311 of file G4RPGReaction.cc.

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}
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cerr
G4DLLIMPORT std::ostream G4cerr

Referenced by NuclearReaction(), and G4RPGTwoCluster::ReactionStage().

◆ GenerateNBodyEventT()

G4double G4RPGReaction::GenerateNBodyEventT ( const G4double  totalEnergy,
const G4bool  constantCrossSection,
std::vector< G4ReactionProduct * > &  list 
)

Definition at line 500 of file G4RPGReaction.cc.

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}

Referenced by G4RPGFragmentation::ReactionStage().

◆ GetFinalStateNucleons()

std::pair< G4int, G4int > G4RPGReaction::GetFinalStateNucleons ( const G4DynamicParticle *  originalTarget,
const G4FastVector< G4ReactionProduct, 256 > &  vec,
const G4int &  vecLen 
)
protected

Definition at line 1011 of file G4RPGReaction.cc.

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 }
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:115

Referenced by G4RPGFragmentation::ReactionStage(), G4RPGTwoBody::ReactionStage(), and G4RPGTwoCluster::ReactionStage().

◆ Isotropic()

G4ThreeVector G4RPGReaction::Isotropic ( const G4double &  pp)
protected

Definition at line 1043 of file G4RPGReaction.cc.

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 }
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42

Referenced by G4RPGFragmentation::ReactionStage(), and G4RPGTwoCluster::ReactionStage().

◆ MomentumCheck()

void G4RPGReaction::MomentumCheck ( const G4ReactionProduct &  modifiedOriginal,
G4ReactionProduct &  currentParticle,
G4ReactionProduct &  targetParticle,
G4FastVector< G4ReactionProduct, 256 > &  vec,
G4int &  vecLen 
)
protected

Definition at line 1054 of file G4RPGReaction.cc.

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 }
G4ReactionProduct::SetTotalEnergy
void SetTotalEnergy(const G4double en)
Definition: G4ReactionProduct.hh:138
G4ReactionProduct::GetTotalMomentum
G4double GetTotalMomentum() const
Definition: G4ReactionProduct.hh:123
G4ReactionProduct::GetMass
G4double GetMass() const
Definition: G4ReactionProduct.hh:147

◆ normal()

G4double G4RPGReaction::normal ( )
protected

Definition at line 685 of file G4RPGReaction.cc.

686{
687 G4double ran = -6.0;
688 for( G4int i=0; i<12; ++i )ran += G4UniformRand();
689 return ran;
690}

Referenced by G4RPGFragmentation::ReactionStage(), G4RPGTwoBody::ReactionStage(), and Rotate().

◆ NuclearReaction()

void G4RPGReaction::NuclearReaction ( G4FastVector< G4ReactionProduct, 4 > &  vec,
G4int &  vecLen,
const G4HadProjectile *  originalIncident,
const G4Nucleus &  aNucleus,
const G4double  theAtomicMass,
const G4double *  massVec 
)

Definition at line 1094 of file G4RPGReaction.cc.

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 }
G4FastVector::Initialize
void Initialize(G4int items)
Definition: G4FastVector.hh:63
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4RPGReaction::GenerateNBodyEvent
G4double GenerateNBodyEvent(const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
Definition: G4RPGReaction.cc:311
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::SetMass
void SetMass(const G4double mas)
Definition: G4ReactionProduct.hh:144

◆ ReactionStage()

G4bool G4RPGReaction::ReactionStage ( const G4HadProjectile *  ,
G4ReactionProduct &  ,
G4bool &  ,
const G4DynamicParticle *  ,
G4ReactionProduct &  ,
G4bool &  ,
const G4Nucleus &  ,
G4ReactionProduct &  ,
G4FastVector< G4ReactionProduct, 256 > &  ,
G4int &  ,
G4bool  ,
G4ReactionProduct &   
)

Definition at line 36 of file G4RPGReaction.cc.

49{
50 G4cout << " G4RPGReactionStage must be overridden in a derived class "
51 << G4endl;
52 return false;
53}
G4cout
G4DLLIMPORT std::ostream G4cout

◆ Rotate()

void G4RPGReaction::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 
)
protected

Definition at line 749 of file G4RPGReaction.cc.

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 }
G4float
float G4float
Definition: G4Types.hh:65
CLHEP::Hep3Vector::cross
Hep3Vector cross(const Hep3Vector &) const
CLHEP::Hep3Vector::dot
double dot(const Hep3Vector &) const
CLHEP::Hep3Vector::rotate
Hep3Vector & rotate(double, const Hep3Vector &)
Definition: ThreeVectorR.cc:28
G4HadProjectile::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4HadProjectile.hh:106
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
G4RPGReaction::Defs1
void Defs1(const G4ReactionProduct &modifiedOriginal, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
Definition: G4RPGReaction.cc:693
G4RPGReaction::normal
G4double normal()
Definition: G4RPGReaction.cc:685
G4INCL::Math::pi
const G4double pi
Definition: G4INCLGlobals.hh:70

Referenced by G4RPGTwoCluster::ReactionStage().

The documentation for this class was generated from the following files: