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

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

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