source: trunk/source/processes/hadronic/stopping/src/G4NeutronCaptureAtRest.cc@ 1036

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

update to geant4.9.2
File size: 12.4 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// G4NeutronCaptureAtRest physics process
27// Larry Felawka (TRIUMF), April 1998
28//---------------------------------------------------------------------
29
30#include "G4NeutronCaptureAtRest.hh"
31#include "G4DynamicParticle.hh"
32#include "G4ParticleTypes.hh"
33#include "Randomize.hh"
34#include <string.h>
35#include <cmath>
36#include <stdio.h>
37
38#define MAX_SECONDARIES 100
39
40// constructor
41
42G4NeutronCaptureAtRest::G4NeutronCaptureAtRest(const G4String& processName,
43 G4ProcessType aType ) :
44 G4VRestProcess (processName, aType), // initialization
45 massProton(G4Proton::Proton()->GetPDGMass()/GeV),
46 massNeutron(G4Neutron::Neutron()->GetPDGMass()/GeV),
47 massElectron(G4Electron::Electron()->GetPDGMass()/GeV),
48 massDeuteron(G4Deuteron::Deuteron()->GetPDGMass()/GeV),
49 massAlpha(G4Alpha::Alpha()->GetPDGMass()/GeV),
50 pdefGamma(G4Gamma::Gamma()),
51 pdefNeutron(G4Neutron::Neutron())
52{
53 if (verboseLevel>0) {
54 G4cout << GetProcessName() << " is created "<< G4endl;
55 }
56 SetProcessSubType(fHadronAtRest);
57 pv = new G4GHEKinematicsVector [MAX_SECONDARIES+1];
58 eve = new G4GHEKinematicsVector [MAX_SECONDARIES];
59 gkin = new G4GHEKinematicsVector [MAX_SECONDARIES];
60
61}
62
63// destructor
64
65G4NeutronCaptureAtRest::~G4NeutronCaptureAtRest()
66{
67 delete [] pv;
68 delete [] eve;
69 delete [] gkin;
70}
71
72
73// methods.............................................................................
74
75G4bool G4NeutronCaptureAtRest::IsApplicable(
76 const G4ParticleDefinition& particle
77 )
78{
79 return ( &particle == pdefNeutron );
80
81}
82
83// Warning - this method may be optimized away if made "inline"
84G4int G4NeutronCaptureAtRest::GetNumberOfSecondaries()
85{
86 return ( ngkine );
87
88}
89
90// Warning - this method may be optimized away if made "inline"
91G4GHEKinematicsVector* G4NeutronCaptureAtRest::GetSecondaryKinematics()
92{
93 return ( &gkin[0] );
94
95}
96
97G4double G4NeutronCaptureAtRest::AtRestGetPhysicalInteractionLength(
98 const G4Track& track,
99 G4ForceCondition* condition
100 )
101{
102 // beggining of tracking
103 ResetNumberOfInteractionLengthLeft();
104
105 // condition is set to "Not Forced"
106 *condition = NotForced;
107
108 // get mean life time
109 currentInteractionLength = GetMeanLifeTime(track, condition);
110
111 if ((currentInteractionLength <0.0) || (verboseLevel>2)){
112 G4cout << "G4NeutronCaptureAtRestProcess::AtRestGetPhysicalInteractionLength ";
113 G4cout << "[ " << GetProcessName() << "]" <<G4endl;
114 track.GetDynamicParticle()->DumpInfo();
115 G4cout << " in Material " << track.GetMaterial()->GetName() <<G4endl;
116 G4cout << "MeanLifeTime = " << currentInteractionLength/ns << "[ns]" <<G4endl;
117 }
118
119 return theNumberOfInteractionLengthLeft * currentInteractionLength;
120
121}
122
123G4VParticleChange* G4NeutronCaptureAtRest::AtRestDoIt(
124 const G4Track& track,
125 const G4Step&
126 )
127//
128// Handles Neutrons at rest; a Neutron can either create secondaries or
129// do nothing (in which case it should be sent back to decay-handling
130// section
131//
132{
133
134// Initialize ParticleChange
135// all members of G4VParticleChange are set to equal to
136// corresponding member in G4Track
137
138 aParticleChange.Initialize(track);
139
140// Store some global quantities that depend on current material and particle
141
142 globalTime = track.GetGlobalTime()/s;
143 G4Material * aMaterial = track.GetMaterial();
144 const G4int numberOfElements = aMaterial->GetNumberOfElements();
145 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
146
147 const G4double* theAtomicNumberDensity = aMaterial->GetAtomicNumDensityVector();
148 G4double normalization = 0;
149 for ( G4int i1=0; i1 < numberOfElements; i1++ )
150 {
151 normalization += theAtomicNumberDensity[i1] ; // change when nucleon specific
152 // probabilities are included.
153 }
154 G4double runningSum= 0.;
155 G4double random = G4UniformRand()*normalization;
156 for ( G4int i2=0; i2 < numberOfElements; i2++ )
157 {
158 runningSum += theAtomicNumberDensity[i2]; // change when nucleon specific
159 // probabilities are included.
160 if (random<=runningSum)
161 {
162 targetCharge = G4double((*theElementVector)[i2]->GetZ());
163 targetAtomicMass = (*theElementVector)[i2]->GetN();
164 }
165 }
166 if (random>runningSum)
167 {
168 targetCharge = G4double((*theElementVector)[numberOfElements-1]->GetZ());
169 targetAtomicMass = (*theElementVector)[numberOfElements-1]->GetN();
170
171 }
172
173 if (verboseLevel>1) {
174 G4cout << "G4NeutronCaptureAtRest::AtRestDoIt is invoked " <<G4endl;
175 }
176
177 G4ParticleMomentum momentum;
178 G4float localtime;
179
180 G4ThreeVector position = track.GetPosition();
181
182 GenerateSecondaries(); // Generate secondaries
183
184 aParticleChange.SetNumberOfSecondaries( ngkine );
185
186 for ( G4int isec = 0; isec < ngkine; isec++ ) {
187 G4DynamicParticle* aNewParticle = new G4DynamicParticle;
188 aNewParticle->SetDefinition( gkin[isec].GetParticleDef() );
189 aNewParticle->SetMomentum( gkin[isec].GetMomentum() * GeV );
190
191 localtime = globalTime + gkin[isec].GetTOF();
192
193 G4Track* aNewTrack = new G4Track( aNewParticle, localtime*s, position );
194 aNewTrack->SetTouchableHandle(track.GetTouchableHandle());
195 aParticleChange.AddSecondary( aNewTrack );
196
197 }
198
199 aParticleChange.ProposeLocalEnergyDeposit( 0.0*GeV );
200
201 aParticleChange.ProposeTrackStatus(fStopAndKill); // Kill the incident Neutron
202
203// clear InteractionLengthLeft
204
205 ResetNumberOfInteractionLengthLeft();
206
207 return &aParticleChange;
208
209}
210
211
212void G4NeutronCaptureAtRest::GenerateSecondaries()
213{
214 static G4int index;
215 static G4int l;
216 static G4int nopt;
217 static G4int i;
218 static G4ParticleDefinition* jnd;
219
220 for (i = 1; i <= MAX_SECONDARIES; ++i) {
221 pv[i].SetZero();
222 }
223
224 ngkine = 0; // number of generated secondary particles
225 ntot = 0;
226 result.SetZero();
227 result.SetMass( massNeutron );
228 result.SetKineticEnergyAndUpdate( 0. );
229 result.SetTOF( 0. );
230 result.SetParticleDef( pdefNeutron );
231
232 NeutronCapture(&nopt);
233
234 // *** CHECK WHETHER THERE ARE NEW PARTICLES GENERATED ***
235 if (ntot != 0 || result.GetParticleDef() != pdefNeutron) {
236 // *** CURRENT PARTICLE IS NOT THE SAME AS IN THE BEGINNING OR/AND ***
237 // *** ONE OR MORE SECONDARIES HAVE BEEN GENERATED ***
238
239 // --- INITIAL PARTICLE TYPE HAS BEEN CHANGED ==> PUT NEW TYPE ON ---
240 // --- THE GEANT TEMPORARY STACK ---
241
242 // --- PUT PARTICLE ON THE STACK ---
243 gkin[0] = result;
244 gkin[0].SetTOF( result.GetTOF() * 5e-11 );
245 ngkine = 1;
246
247 // --- ALL QUANTITIES ARE TAKEN FROM THE GHEISHA STACK WHERE THE ---
248 // --- CONVENTION IS THE FOLLOWING ---
249
250 // --- ONE OR MORE SECONDARIES HAVE BEEN GENERATED ---
251 for (l = 1; l <= ntot; ++l) {
252 index = l - 1;
253 jnd = eve[index].GetParticleDef();
254
255 // --- ADD PARTICLE TO THE STACK IF STACK NOT YET FULL ---
256 if (ngkine < MAX_SECONDARIES) {
257 gkin[ngkine] = eve[index];
258 gkin[ngkine].SetTOF( eve[index].GetTOF() * 5e-11 );
259 ++ngkine;
260 }
261 }
262 }
263 else {
264 // --- NO SECONDARIES GENERATED AND PARTICLE IS STILL THE SAME ---
265 // --- ==> COPY EVERYTHING BACK IN THE CURRENT GEANT STACK ---
266 ngkine = 0;
267 ntot = 0;
268 globalTime += result.GetTOF() * G4float(5e-11);
269 }
270
271 // --- LIMIT THE VALUE OF NGKINE IN CASE OF OVERFLOW ---
272 ngkine = G4int(std::min(ngkine,G4int(MAX_SECONDARIES)));
273
274} // GenerateSecondaries
275
276
277void G4NeutronCaptureAtRest::Normal(G4float *ran)
278{
279 static G4int i;
280
281 // *** NVE 14-APR-1988 CERN GENEVA ***
282 // ORIGIN : H.FESEFELDT (27-OCT-1983)
283
284 *ran = G4float(-6.);
285 for (i = 1; i <= 12; ++i) {
286 *ran += G4UniformRand();
287 }
288
289} // Normal
290
291
292void G4NeutronCaptureAtRest::NeutronCapture(G4int *nopt)
293{
294 static G4int nt;
295 static G4float xp, pcm;
296 static G4float ran;
297
298 // *** ROUTINE FOR CAPTURE OF NEUTRAL BARYONS ***
299 // *** NVE 04-MAR-1988 CERN GENEVA ***
300 // ORIGIN : H.FESEFELDT (02-DEC-1986)
301
302 *nopt = 1;
303 pv[1] = result;
304 pv[2].SetZero();
305 pv[2].SetMass( AtomAs(targetAtomicMass, targetCharge) );
306 pv[2].SetMomentumAndUpdate( 0., 0., 0. );
307 pv[2].SetTOF( result.GetTOF() );
308 pv[2].SetParticleDef( NULL );
309 pv[MAX_SECONDARIES].Add( pv[1], pv[2] );
310 pv[MAX_SECONDARIES].SetMomentum( -pv[MAX_SECONDARIES].GetMomentum().x(), -pv[MAX_SECONDARIES].GetMomentum().y(), -pv[MAX_SECONDARIES].GetMomentum().z() );
311 pv[MAX_SECONDARIES].SetParticleDef( NULL );
312 Normal(&ran);
313 pcm = ran * G4float(.001) + G4float(.0065);
314 ran = G4UniformRand();
315 result.SetTOF( result.GetTOF() - std::log(ran) * G4float(480.) );
316 pv[3].SetZero();
317 pv[3].SetMass( 0. );
318 pv[3].SetKineticEnergyAndUpdate( pcm );
319 pv[3].SetTOF( result.GetTOF() );
320 pv[3].SetParticleDef( pdefGamma );
321 pv[3].Lor( pv[3], pv[MAX_SECONDARIES] );
322 nt = 3;
323 xp = G4float(.008) - pcm;
324 if (xp >= G4float(0.)) {
325 nt = 4;
326 pv[4].SetZero();
327 pv[4].SetMass( 0. );
328 pv[4].SetKineticEnergyAndUpdate( xp );
329 pv[4].SetTOF( result.GetTOF() );
330 pv[4].SetParticleDef( pdefGamma );
331 pv[4].Lor( pv[4], pv[MAX_SECONDARIES] );
332 }
333 result = pv[3];
334 if (nt == 4) {
335 if (ntot < MAX_SECONDARIES-1) {
336 eve[ntot++] = pv[4];
337 }
338 }
339
340} // NeutronCapture
341
342
343G4double G4NeutronCaptureAtRest::AtomAs(G4float a, G4float z)
344{
345 G4float ret_val;
346 G4double d__1, d__2;
347
348 static G4double aa;
349 static G4int ia, iz;
350 static G4double zz;
351 static G4float rma, rmd;
352 static G4int ipp;
353 static G4float rmn, rmp;
354 static G4int izz;
355 static G4float rmel;
356 static G4double mass;
357
358 // *** DETERMINATION OF THE ATOMIC MASS ***
359 // *** NVE 19-MAY-1988 CERN GENEVA ***
360 // ORIGIN : H.FESEFELDT (02-DEC-1986)
361
362 // --- GET ATOMIC (= ELECTRONS INCL.) MASSES (IN MEV) FROM RMASS ARRAY ---
363 // --- ELECTRON ---
364 rmel = massElectron * G4float(1e3);
365 // --- PROTON ---
366 rmp = massProton * G4float(1e3);
367 // --- NEUTRON ---
368 rmn = massNeutron * G4float(1e3);
369 // --- DEUTERON ---
370 rmd = massDeuteron * G4float(1e3) + rmel;
371 // --- ALPHA ---
372 rma = massAlpha * G4float(1e3) + rmel * G4float(2.);
373
374 ret_val = G4float(0.);
375 aa = a * 1.;
376 zz = z * 1.;
377 ia = G4int(a + G4float(.5));
378 if (ia < 1) {
379 return ret_val;
380 }
381 iz = G4int(z + G4float(.5));
382 if (iz < 0 || iz > ia) {
383 return ret_val;
384 }
385 mass = 0.;
386 if (ia == 1) {
387 if (iz == 0) {
388 mass = rmn;
389 }
390 else if (iz == 1) {
391 mass = rmp + rmel;
392 }
393 }
394 else if (ia == 2 && iz == 1) {
395 mass = rmd;
396 }
397 else if (ia == 4 && iz == 2) {
398 mass = rma;
399 }
400 else if ( (ia == 2 && iz != 1) || ia == 3 || (ia == 4 && iz != 2) || ia > 4) {
401 d__1 = aa / G4float(2.) - zz;
402 d__2 = zz;
403 mass = (aa - zz) * rmn + zz * rmp + zz * rmel - aa * G4float(15.67) +
404 std::pow(aa, .6666667) * G4float(17.23) + d__1 * d__1 * G4float(93.15) / aa +
405 d__2 * d__2 * G4float(.6984523) / std::pow(aa, .3333333);
406 ipp = (ia - iz) % 2;
407 izz = iz % 2;
408 if (ipp == izz) {
409 mass += (ipp + izz - 1) * G4float(12.) * std::pow(aa, -.5);
410 }
411 }
412 ret_val = mass * G4float(.001);
413 return ret_val;
414
415} // AtomAs
Note: See TracBrowser for help on using the repository browser.