source: trunk/source/processes/hadronic/models/rpg/src/G4RPGProtonInelastic.cc@ 819

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

import all except CVS
File size: 15.9 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: G4RPGProtonInelastic.cc,v 1.1 2007/07/18 21:04:20 dennis Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-02 $
28//
29
30#include "G4RPGProtonInelastic.hh"
31#include "Randomize.hh"
32
33G4HadFinalState*
34G4RPGProtonInelastic::ApplyYourself( const G4HadProjectile &aTrack,
35 G4Nucleus &targetNucleus )
36{
37 theParticleChange.Clear();
38 const G4HadProjectile *originalIncident = &aTrack;
39 if (originalIncident->GetKineticEnergy()<= 0.1*MeV)
40 {
41 theParticleChange.SetStatusChange(isAlive);
42 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
43 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
44 return &theParticleChange;
45 }
46
47 //
48 // create the target particle
49 //
50 G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
51 if( verboseLevel > 1 )
52 {
53 const G4Material *targetMaterial = aTrack.GetMaterial();
54 G4cout << "G4RPGProtonInelastic::ApplyYourself called" << G4endl;
55 G4cout << "kinetic energy = "
56 << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
57 G4cout << "target material = "
58 << targetMaterial->GetName() << ", ";
59 G4cout << "target particle = "
60 << originalTarget->GetDefinition()->GetParticleName()
61 << G4endl;
62 }
63
64 if( originalIncident->GetKineticEnergy()/GeV < 0.01+2.*G4UniformRand()/9. )
65 {
66 SlowProton( originalIncident, targetNucleus );
67 delete originalTarget;
68 return &theParticleChange;
69 }
70
71 //
72 // Fermi motion and evaporation
73 // As of Geant3, the Fermi energy calculation had not been Done
74 //
75 G4double ek = originalIncident->GetKineticEnergy()/MeV;
76 G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
77 G4ReactionProduct modifiedOriginal;
78 modifiedOriginal = *originalIncident;
79
80 G4double tkin = targetNucleus.Cinema( ek );
81 ek += tkin;
82 modifiedOriginal.SetKineticEnergy( ek*MeV );
83 G4double et = ek + amas;
84 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
85 G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
86 if( pp > 0.0 )
87 {
88 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
89 modifiedOriginal.SetMomentum( momentum * (p/pp) );
90 }
91 //
92 // calculate black track energies
93 //
94 tkin = targetNucleus.EvaporationEffects( ek );
95 ek -= tkin;
96 modifiedOriginal.SetKineticEnergy( ek*MeV );
97 et = ek + amas;
98 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
99 pp = modifiedOriginal.GetMomentum().mag()/MeV;
100 if( pp > 0.0 )
101 {
102 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
103 modifiedOriginal.SetMomentum( momentum * (p/pp) );
104 }
105 const G4double cutOff = 0.1;
106 if( modifiedOriginal.GetKineticEnergy()/MeV <= cutOff )
107 {
108 SlowProton( originalIncident, targetNucleus );
109 delete originalTarget;
110 return &theParticleChange;
111 }
112
113 G4ReactionProduct currentParticle = modifiedOriginal;
114 G4ReactionProduct targetParticle;
115 targetParticle = *originalTarget;
116 currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
117 targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
118 G4bool incidentHasChanged = false;
119 G4bool targetHasChanged = false;
120 G4bool quasiElastic = false;
121 G4FastVector<G4ReactionProduct,256> vec; // vec will contain the sec. particles
122 G4int vecLen = 0;
123 vec.Initialize( 0 );
124
125 Cascade( vec, vecLen,
126 originalIncident, currentParticle, targetParticle,
127 incidentHasChanged, targetHasChanged, quasiElastic );
128
129 CalculateMomenta( vec, vecLen,
130 originalIncident, originalTarget, modifiedOriginal,
131 targetNucleus, currentParticle, targetParticle,
132 incidentHasChanged, targetHasChanged, quasiElastic );
133
134 SetUpChange( vec, vecLen,
135 currentParticle, targetParticle,
136 incidentHasChanged );
137
138 delete originalTarget;
139 return &theParticleChange;
140}
141
142void
143G4RPGProtonInelastic::SlowProton(const G4HadProjectile *originalIncident,
144 G4Nucleus &targetNucleus )
145{
146 const G4double A = targetNucleus.GetN(); // atomic weight
147 const G4double Z = targetNucleus.GetZ(); // atomic number
148// G4double currentKinetic = originalIncident->GetKineticEnergy()/MeV;
149 //
150 // calculate Q-value of reactions
151 //
152 G4double theAtomicMass = targetNucleus.AtomicMass( A, Z );
153 G4double massVec[9];
154 massVec[0] = targetNucleus.AtomicMass( A+1.0, Z+1.0 );
155 massVec[1] = 0.;
156 if (A > Z+1.0)
157 massVec[1] = targetNucleus.AtomicMass( A , Z+1.0 );
158 massVec[2] = theAtomicMass;
159 massVec[3] = 0.;
160 if (A > 1.0 && A-1.0 > Z)
161 massVec[3] = targetNucleus.AtomicMass( A-1.0, Z );
162 massVec[4] = 0.;
163 if (A > 2.0 && A-2.0 > Z)
164 massVec[4] = targetNucleus.AtomicMass( A-2.0, Z );
165 massVec[5] = 0.;
166 if (A > 3.0 && Z > 1.0 && A-3.0 > Z-1.0)
167 massVec[5] = targetNucleus.AtomicMass( A-3.0, Z-1.0 );
168 massVec[6] = 0.;
169 if (A > 1.0 && A-1.0 > Z+1.0)
170 massVec[6] = targetNucleus.AtomicMass( A-1.0, Z+1.0 );
171 massVec[7] = massVec[3];
172 massVec[8] = 0.;
173 if (A > 1.0 && Z > 1.0)
174 massVec[8] = targetNucleus.AtomicMass( A-1.0, Z-1.0 );
175
176 G4FastVector<G4ReactionProduct,4> vec; // vec will contain the secondary particles
177 G4int vecLen = 0;
178 vec.Initialize( 0 );
179
180 twoBody.NuclearReaction( vec, vecLen, originalIncident,
181 targetNucleus, theAtomicMass, massVec );
182
183 theParticleChange.SetStatusChange( stopAndKill );
184 theParticleChange.SetEnergyChange( 0.0 );
185
186 G4DynamicParticle *pd;
187 for( G4int i=0; i<vecLen; ++i )
188 {
189 pd = new G4DynamicParticle();
190 pd->SetDefinition( vec[i]->GetDefinition() );
191 pd->SetMomentum( vec[i]->GetMomentum() );
192 theParticleChange.AddSecondary( pd );
193 delete vec[i];
194 }
195}
196
197void G4RPGProtonInelastic::Cascade(
198 G4FastVector<G4ReactionProduct,256> &vec,
199 G4int &vecLen,
200 const G4HadProjectile *originalIncident,
201 G4ReactionProduct &currentParticle,
202 G4ReactionProduct &targetParticle,
203 G4bool &incidentHasChanged,
204 G4bool &targetHasChanged,
205 G4bool &quasiElastic )
206{
207 // Derived from H. Fesefeldt's original FORTRAN code CASP
208 //
209 // Proton undergoes interaction with nucleon within a nucleus. Check if
210 // it is energetically possible to produce pions/kaons. In not, assume
211 // nuclear excitation occurs and input particle is degraded in energy. No
212 // other particles are produced.
213 // If reaction is possible, find the correct number of
214 // pions/protons/neutrons produced using an interpolation to multiplicity
215 // data. Replace some pions or protons/neutrons by kaons or strange
216 // baryons according to the average multiplicity per Inelastic reaction.
217 //
218 // The center of mass energy is based on the initial energy, before
219 // Fermi motion and evaporation effects are taken into account
220 //
221 const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
222 const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
223 const G4double targetMass = targetParticle.GetMass()/MeV;
224 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
225 targetMass*targetMass +
226 2.0*targetMass*etOriginal );
227 G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
228 if( availableEnergy <= G4PionPlus::PionPlus()->GetPDGMass()/MeV )
229 { // not energetically possible to produce pion(s)
230 quasiElastic = true;
231 return;
232 }
233 static G4bool first = true;
234 const G4int numMul = 1200;
235 const G4int numSec = 60;
236 static G4double protmul[numMul], protnorm[numSec]; // proton constants
237 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
238 // np = number of pi+, nm = number of pi-, nz = number of pi0
239 G4int counter, nt=0, np=0, nm=0, nz=0;
240 const G4double c = 1.25;
241 const G4double b[] = { 0.70, 0.35 };
242 if( first ) // compute normalization constants, this will only be Done once
243 {
244 first = false;
245 G4int i;
246 for( i=0; i<numMul; ++i )protmul[i] = 0.0;
247 for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
248 counter = -1;
249 for( np=0; np<(numSec/3); ++np )
250 {
251 for( nm=std::max(0,np-2); nm<=np; ++nm )
252 {
253 for( nz=0; nz<numSec/3; ++nz )
254 {
255 if( ++counter < numMul )
256 {
257 nt = np+nm+nz;
258 if( nt>0 && nt<=numSec )
259 {
260 protmul[counter] = Pmltpc(np,nm,nz,nt,b[0],c) /
261 (Factorial(2-np+nm)*Factorial(np-nm) );
262 protnorm[nt-1] += protmul[counter];
263 }
264 }
265 }
266 }
267 }
268 for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
269 for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
270 counter = -1;
271 for( np=0; np<numSec/3; ++np )
272 {
273 for( nm=std::max(0,np-1); nm<=(np+1); ++nm )
274 {
275 for( nz=0; nz<numSec/3; ++nz )
276 {
277 if( ++counter < numMul )
278 {
279 nt = np+nm+nz;
280 if( nt>0 && nt<=numSec )
281 {
282 neutmul[counter] = Pmltpc(np,nm,nz,nt,b[1],c) /
283 (Factorial(1-np+nm)*Factorial(1+np-nm) );
284 neutnorm[nt-1] += neutmul[counter];
285 }
286 }
287 }
288 }
289 }
290 for( i=0; i<numSec; ++i )
291 {
292 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
293 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
294 }
295 } // end of initialization
296
297 const G4double expxu = 82.; // upper bound for arg. of exp
298 const G4double expxl = -expxu; // lower bound for arg. of exp
299 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
300 G4ParticleDefinition *aProton = G4Proton::Proton();
301 G4int ieab = static_cast<G4int>(availableEnergy*5.0/GeV);
302 const G4double supp[] = {0.,0.4,0.55,0.65,0.75,0.82,0.86,0.90,0.94,0.98};
303 G4double test, w0, wp, wt, wm;
304 if( (availableEnergy < 2.0*GeV) && (G4UniformRand() >= supp[ieab]) )
305 {
306 // suppress high multiplicity events at low momentum
307 // only one pion will be produced
308
309 np = nm = nz = 0;
310 if( targetParticle.GetDefinition() == aProton )
311 {
312 test = std::exp( std::min( expxu, std::max(
313 expxl, -(1.0+b[0])*(1.0+b[0])/(2.0*c*c) ) ) );
314 w0 = test/2.0;
315 wp = test;
316 if( G4UniformRand() < w0/(w0+wp) )
317 nz = 1;
318 else
319 np = 1;
320 }
321 else // target is a neutron
322 {
323 test = std::exp( std::min( expxu, std::max(
324 expxl, -(1.0+b[1])*(1.0+b[1])/(2.0*c*c) ) ) );
325 w0 = test;
326 wp = test/2.0;
327 test = std::exp( std::min( expxu, std::max(
328 expxl, -(-1.0+b[1])*(-1.0+b[1])/(2.0*c*c) ) ) );
329 wm = test/2.0;
330 wt = w0+wp+wm;
331 wp += w0;
332 G4double ran = G4UniformRand();
333 if( ran < w0/wt )
334 nz = 1;
335 else if( ran < wp/wt )
336 np = 1;
337 else
338 nm = 1;
339 }
340 }
341 else // (availableEnergy >= 2.0*GeV) || (random number < supp[ieab])
342 {
343 G4double n, anpn;
344 GetNormalizationConstant( availableEnergy, n, anpn );
345 G4double ran = G4UniformRand();
346 G4double dum, excs = 0.0;
347 if( targetParticle.GetDefinition() == aProton )
348 {
349 counter = -1;
350 for( np=0; np<numSec/3 && ran>=excs; ++np )
351 {
352 for( nm=std::max(0,np-2); nm<=np && ran>=excs; ++nm )
353 {
354 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
355 {
356 if( ++counter < numMul )
357 {
358 nt = np+nm+nz;
359 if( nt>0 && nt<=numSec )
360 {
361 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
362 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
363 if( std::fabs(dum) < 1.0 )
364 {
365 if( test >= 1.0e-10 )excs += dum*test;
366 } else {
367 excs += dum*test;
368 }
369 }
370 }
371 }
372 }
373 }
374 if( ran >= excs ) // 3 previous loops continued to the end
375 {
376 quasiElastic = true;
377 return;
378 }
379 }
380 else // target must be a neutron
381 {
382 counter = -1;
383 for( np=0; np<numSec/3 && ran>=excs; ++np )
384 {
385 for( nm=std::max(0,np-1); nm<=(np+1) && ran>=excs; ++nm )
386 {
387 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
388 {
389 if( ++counter < numMul )
390 {
391 nt = np+nm+nz;
392 if( nt>0 && nt<=numSec )
393 {
394 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
395 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
396 if( std::fabs(dum) < 1.0 )
397 {
398 if( test >= 1.0e-10 )excs += dum*test;
399 } else {
400 excs += dum*test;
401 }
402 }
403 }
404 }
405 }
406 }
407 if( ran >= excs ) // 3 previous loops continued to the end
408 {
409 quasiElastic = true;
410 return;
411 }
412 }
413 np--; nm--; nz--;
414 }
415 if( targetParticle.GetDefinition() == aProton )
416 {
417 switch( np-nm )
418 {
419 case 1:
420 if( G4UniformRand() < 0.5 )
421 {
422 targetParticle.SetDefinitionAndUpdateE( aNeutron );
423 targetHasChanged = true;
424 } else {
425 currentParticle.SetDefinitionAndUpdateE( aNeutron );
426 incidentHasChanged = true;
427 }
428 break;
429 case 2:
430 currentParticle.SetDefinitionAndUpdateE( aNeutron );
431 targetParticle.SetDefinitionAndUpdateE( aNeutron );
432 incidentHasChanged = true;
433 targetHasChanged = true;
434 break;
435 default:
436 break;
437 }
438 }
439 else // target is a neutron
440 {
441 switch( np-nm )
442 {
443 case 0:
444 if( G4UniformRand() < 0.333333 )
445 {
446 currentParticle.SetDefinitionAndUpdateE( aNeutron );
447 targetParticle.SetDefinitionAndUpdateE( aProton );
448 incidentHasChanged = true;
449 targetHasChanged = true;
450 }
451 break;
452 case 1:
453 currentParticle.SetDefinitionAndUpdateE( aNeutron );
454 incidentHasChanged = true;
455 break;
456 default:
457 targetParticle.SetDefinitionAndUpdateE( aProton );
458 targetHasChanged = true;
459 break;
460 }
461 }
462 SetUpPions( np, nm, nz, vec, vecLen );
463 return;
464 }
465
466 /* end of file */
467
Note: See TracBrowser for help on using the repository browser.