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

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

import all except CVS
File size: 12.3 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
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.