source: trunk/source/processes/hadronic/models/de_excitation/ablation/src/G4WilsonAblationModel.cc @ 819

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

import all except CVS
File size: 16.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// *                                                                  *
21// * Parts of this code which have been  developed by QinetiQ Ltd     *
22// * under contract to the European Space Agency (ESA) are the        *
23// * intellectual property of ESA. Rights to use, copy, modify and    *
24// * redistribute this software for general public use are granted    *
25// * in compliance with any licensing, distribution and development   *
26// * policy adopted by the Geant4 Collaboration. This code has been   *
27// * written by QinetiQ Ltd for the European Space Agency, under ESA  *
28// * contract 17191/03/NL/LvH (Aurora Programme).                     *
29// *                                                                  *
30// * By using,  copying,  modifying or  distributing the software (or *
31// * any work based  on the software)  you  agree  to acknowledge its *
32// * use  in  resulting  scientific  publications,  and indicate your *
33// * acceptance of all terms of the Geant4 Software license.          *
34// ********************************************************************
35//
36// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37//
38// MODULE:              G4WilsonAblationModel.cc
39//
40// Version:             B.1
41// Date:                15/04/04
42// Author:              P R Truscott
43// Organisation:        QinetiQ Ltd, UK
44// Customer:            ESA/ESTEC, NOORDWIJK
45// Contract:            17191/03/NL/LvH
46//
47// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48//
49// CHANGE HISTORY
50// --------------
51//
52// 6 October 2003, P R Truscott, QinetiQ Ltd, UK
53// Created.
54//
55// 15 March 2004, P R Truscott, QinetiQ Ltd, UK
56// Beta release
57//
58// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
59////////////////////////////////////////////////////////////////////////////////
60//
61#include "G4WilsonAblationModel.hh"
62#include "Randomize.hh"
63#include "G4ParticleTable.hh"
64#include "G4IonTable.hh"
65#include "G4Alpha.hh"
66#include "G4He3.hh"
67#include "G4Triton.hh"
68#include "G4Deuteron.hh"
69#include "G4Proton.hh"
70#include "G4Neutron.hh"
71#include "G4AlphaEvaporationChannel.hh"
72#include "G4He3EvaporationChannel.hh"
73#include "G4TritonEvaporationChannel.hh"
74#include "G4DeuteronEvaporationChannel.hh"
75#include "G4ProtonEvaporationChannel.hh"
76#include "G4NeutronEvaporationChannel.hh"
77#include "G4LorentzVector.hh"
78#include "G4VEvaporationChannel.hh"
79
80#include <iomanip>
81#include <numeric>
82////////////////////////////////////////////////////////////////////////////////
83//
84G4WilsonAblationModel::G4WilsonAblationModel()
85{
86//
87//
88// Send message to stdout to advise that the G4Abrasion model is being used.
89//
90  PrintWelcomeMessage();
91//
92//
93// Set the default verbose level to 0 - no output.
94//
95  verboseLevel = 0; 
96//
97//
98// Set the binding energy per nucleon .... did I mention that this is a crude
99// model for nuclear de-excitation?
100//
101  B = 10.0 * MeV;
102//
103//
104// It is possuble to switch off secondary particle production (other than the
105// final nuclear fragment).  The default is on.
106//
107  produceSecondaries = true;
108//
109//
110// Now we need to define the decay modes.  We're using the G4Evaporation model
111// to help determine the kinematics of the decay.
112//
113  nFragTypes  = 6;
114  fragType[0] = G4Alpha::Alpha();
115  fragType[1] = G4He3::He3();
116  fragType[2] = G4Triton::Triton();
117  fragType[3] = G4Deuteron::Deuteron();
118  fragType[4] = G4Proton::Proton();
119  fragType[5] = G4Neutron::Neutron();
120//
121//
122// Set verboseLevel default to no output.
123//
124  verboseLevel = 0;
125}
126////////////////////////////////////////////////////////////////////////////////
127//
128G4WilsonAblationModel::~G4WilsonAblationModel()
129{;}
130////////////////////////////////////////////////////////////////////////////////
131//
132G4FragmentVector *G4WilsonAblationModel::BreakItUp
133  (const G4Fragment &theNucleus)
134{
135//
136//
137// Initilise the pointer to the G4FragmentVector used to return the information
138// about the breakup.
139//
140  fragmentVector = new G4FragmentVector;
141  fragmentVector->clear();
142//
143//
144// Get the A, Z and excitation of the nucleus.
145//
146  G4int A     = (G4int) theNucleus.GetA();
147  G4int Z     = (G4int) theNucleus.GetZ();
148  G4double ex = theNucleus.GetExcitationEnergy();
149  if (verboseLevel >= 2)
150  {
151    G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
152           <<"oooooooooooooooooooooooooooooooooooooooo"
153           <<G4endl;
154    G4cout.precision(6);
155    G4cout <<"IN G4WilsonAblationModel" <<G4endl;
156    G4cout <<"Initial prefragment A=" <<A
157           <<", Z=" <<Z
158           <<", excitation energy = " <<ex/MeV <<" MeV"
159           <<G4endl; 
160  }
161//
162//
163// Check that there is a nucleus to speak of.  It's possible there isn't one
164// or its just a proton or neutron.  In either case, the excitation energy
165// (from the Lorentz vector) is not used.
166//
167  if (A == 0)
168  {
169    if (verboseLevel >= 2)
170    {
171      G4cout <<"No nucleus to decay" <<G4endl;
172      G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
173             <<"oooooooooooooooooooooooooooooooooooooooo"
174             <<G4endl;
175    }
176    return fragmentVector;
177  }
178  else if (A == 1)
179  {
180    G4LorentzVector lorentzVector = theNucleus.GetMomentum();
181    lorentzVector.setE(lorentzVector.e()-ex+10.0*eV);
182    if (Z == 0)
183    {
184      G4Fragment *fragment = new G4Fragment(lorentzVector,G4Neutron::Neutron());
185      fragmentVector->push_back(fragment);
186    }
187    else
188    {
189      G4Fragment *fragment = new G4Fragment(lorentzVector,G4Proton::Proton());
190      fragmentVector->push_back(fragment);
191    }
192    if (verboseLevel >= 2)
193    {
194      G4cout <<"Final fragment is in fact only a nucleon) :" <<G4endl;
195      G4cout <<(*fragmentVector)[0] <<G4endl;
196      G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
197             <<"oooooooooooooooooooooooooooooooooooooooo"
198             <<G4endl;
199    }
200    return fragmentVector;
201  }
202//
203//
204// Then the number of nucleons ablated (either as nucleons or light nuclear
205// fragments) is based on a simple argument for the binding energy per nucleon.
206//
207  G4int DAabl = (G4int) (ex / B);
208  if (DAabl > A) DAabl = A;
209  if (verboseLevel >= 2)
210    G4cout <<"Number of nucleons ejected = " <<DAabl <<G4endl;
211
212//
213//
214// Determine the nuclear fragment from the ablation process by sampling the
215// Rudstam equation.
216//
217  G4int AF = A - DAabl;
218  G4int ZF = 0;
219  if (AF > 0)
220  {
221    G4double AFd = static_cast<G4double>(AF);
222    G4double R = 11.8 / std::pow(AFd, 0.45);
223    G4int minZ = Z - DAabl;
224    if (minZ <= 0) minZ = 1;
225//
226//
227// Here we define an integral probability distribution based on the Rudstam
228// equation assuming a constant AF.
229//   
230    G4double sig[100];
231    G4double sum = 0.0;
232    for (G4int ii=minZ; ii<= Z; ii++)
233    {
234      sum   += std::exp(-R*std::pow(std::abs(ii - 0.486*AFd + 3.8E-04*AFd*AFd),1.5));
235      sig[ii] = sum;
236    }
237//
238//
239// Now sample that distribution to determine a value for ZF.
240//
241    G4double xi  = G4UniformRand();
242    G4int iz     = minZ;
243    G4bool found = false;
244    while (iz <= Z && !found)
245    {
246      found = (xi <= sig[iz]/sum);
247      if (!found) iz++;
248    }
249    if (iz > Z)
250      ZF = Z;
251    else
252      ZF = iz;
253  }
254  G4int DZabl = Z - ZF;
255  if (verboseLevel >= 2)
256    G4cout <<"Final fragment      A=" <<AF
257           <<", Z=" <<ZF
258           <<G4endl;
259//
260//
261// Now determine the nucleons or nuclei which have bee ablated.  The preference
262// is for the production of alphas, then other nuclei in order of decreasing
263// binding energy. The energies assigned to the products of the decay are
264// provisional for the moment (the 10eV is just to avoid errors with negative
265// excitation energies due to rounding).
266//
267  G4double totalEpost = 0.0;
268  evapType.clear();
269  for (G4int ift=0; ift<nFragTypes; ift++)
270  {
271    G4ParticleDefinition *type = fragType[ift];
272    G4double n  = std::floor((G4double) DAabl / type->GetBaryonNumber() + 1.0E-10);
273    G4double n1 = 1.0E+10;
274    if (fragType[ift]->GetPDGCharge() > 0.0)
275      n1 = std::floor((G4double) DZabl / type->GetPDGCharge() + 1.0E-10);
276    if (n > n1) n = n1;
277    if (n > 0.0)
278    {
279      G4double mass = type->GetPDGMass();
280      for (G4int j=0; j<(G4int) n; j++)
281      {
282        totalEpost += mass;
283        evapType.push_back(type);
284      }
285      DAabl -= (G4int) (n * type->GetBaryonNumber() + 1.0E-10);
286      DZabl -= (G4int) (n * type->GetPDGCharge() + 1.0E-10);
287      if (verboseLevel >= 2)
288        G4cout <<"Particle type: " <<std::setw(10) <<type->GetParticleName()
289               <<", number of particles emitted = " <<n
290               <<G4endl;
291    }
292  }
293//
294//
295// Determine the properties of the final nuclear fragment.
296//
297  G4double massFinalFrag = 0.0;
298  if (AF > 0.0)
299    massFinalFrag = G4ParticleTable::GetParticleTable()->GetIonTable()->
300      GetIonMass(ZF,AF);
301  totalEpost   += massFinalFrag;
302//
303//
304// Add the total energy from the fragment.  Note that the fragment is assumed
305// to be de-excited and does not undergo photo-evaporation .... I did mention
306// this is a bit of a crude model?
307//
308  G4double massPreFrag      = theNucleus.GetGroundStateMass();
309  G4double totalEpre        = massPreFrag + ex;
310  G4double excess           = totalEpre - totalEpost;
311//  G4Fragment *resultNucleus(theNucleus);
312  G4Fragment *resultNucleus = new G4Fragment(A, Z, theNucleus.GetMomentum());
313  G4ThreeVector boost(0.0,0.0,0.0);
314  G4int nEvap = 0;
315  if (produceSecondaries && evapType.size()>0)
316  {
317    if (excess > 0.0)
318    {
319      SelectSecondariesByEvaporation (resultNucleus);
320      nEvap = fragmentVector->size();
321      boost = resultNucleus->GetMomentum().findBoostToCM();
322      if (evapType.size() > 0)
323        SelectSecondariesByDefault (boost);
324    }
325    else
326      SelectSecondariesByDefault(G4ThreeVector(0.0,0.0,0.0));
327  }
328  if (AF > 0)
329  {
330    G4double mass = G4ParticleTable::GetParticleTable()->GetIonTable()->
331      GetIonMass(ZF,AF);
332    G4double e    = mass + 10.0*eV;
333    G4double p    = std::sqrt(e*e-mass*mass);
334    G4ThreeVector direction(0.0,0.0,1.0);
335    G4LorentzVector lorentzVector = G4LorentzVector(direction*p, e);
336    lorentzVector.boost(-boost);
337    *resultNucleus = G4Fragment(AF, ZF, lorentzVector);
338    fragmentVector->push_back(resultNucleus);
339  }
340//
341//
342// Provide verbose output on the ablation products if requested.
343//
344  if (verboseLevel >= 2)
345  {
346    if (nEvap > 0)
347    {
348      G4cout <<"----------------------" <<G4endl;
349      G4cout <<"Evaporated particles :" <<G4endl;
350      G4cout <<"----------------------" <<G4endl;
351    }
352    G4int ie = 0;
353    G4FragmentVector::iterator iter;
354    for (iter = fragmentVector->begin(); iter != fragmentVector->end(); ++iter)
355    {
356      if (ie == nEvap)
357      {
358        G4cout <<*iter  <<G4endl;
359        G4cout <<"---------------------------------" <<G4endl;
360        G4cout <<"Particles from default emission :" <<G4endl;
361        G4cout <<"---------------------------------" <<G4endl;
362      }
363      G4cout <<*iter <<G4endl;
364    }
365    G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
366           <<"oooooooooooooooooooooooooooooooooooooooo"
367           <<G4endl;
368  }
369
370  return fragmentVector;   
371}
372////////////////////////////////////////////////////////////////////////////////
373//
374void G4WilsonAblationModel::SelectSecondariesByEvaporation
375  (G4Fragment *intermediateNucleus)
376{
377  G4bool evaporate = true;
378  while (evaporate && evapType.size() != 0)
379  {
380//
381//
382// Here's the cheaky bit.  We're hijacking the G4Evaporation model, in order to
383// more accurately sample to kinematics, but the species of the nuclear
384// fragments will be the ones of our choosing as above.
385//
386    std::vector <G4VEvaporationChannel*>  theChannels;
387    theChannels.clear();
388    VectorOfFragmentTypes::iterator iter;
389    std::vector <VectorOfFragmentTypes::iterator> iters;
390    iters.clear();
391    iter = std::find(evapType.begin(), evapType.end(), G4Alpha::Alpha());
392    if (iter != evapType.end())
393    {
394      theChannels.push_back(new G4AlphaEvaporationChannel);
395      iters.push_back(iter);
396    }
397    iter = std::find(evapType.begin(), evapType.end(), G4He3::He3());
398    if (iter != evapType.end())
399    {
400      theChannels.push_back(new G4He3EvaporationChannel);
401      iters.push_back(iter);
402    }
403    iter = std::find(evapType.begin(), evapType.end(), G4Triton::Triton());
404    if (iter != evapType.end())
405    {
406      theChannels.push_back(new G4TritonEvaporationChannel);
407      iters.push_back(iter);
408    }
409    iter = std::find(evapType.begin(), evapType.end(), G4Deuteron::Deuteron());
410    if (iter != evapType.end())
411    {
412      theChannels.push_back(new G4DeuteronEvaporationChannel);
413      iters.push_back(iter);
414    }
415    iter = std::find(evapType.begin(), evapType.end(), G4Proton::Proton());
416    if (iter != evapType.end())
417    {
418      theChannels.push_back(new G4ProtonEvaporationChannel);
419      iters.push_back(iter);
420    }
421    iter = std::find(evapType.begin(), evapType.end(), G4Neutron::Neutron());
422    if (iter != evapType.end())
423    {
424      theChannels.push_back(new G4NeutronEvaporationChannel);
425      iters.push_back(iter);
426    }
427    G4int nChannels = theChannels.size();
428
429    std::vector<G4VEvaporationChannel*>::iterator iterEv;
430    for (iterEv=theChannels.begin(); iterEv!=theChannels.end(); iterEv++)
431      (*iterEv)->Initialize(*intermediateNucleus);
432    G4double totalProb = std::accumulate(theChannels.begin(),
433      theChannels.end(), 0.0, SumProbabilities());
434    if (totalProb > 0.0)
435    {
436//
437//
438// The emission probability for at least one of the evaporation channels is
439// positive, therefore work out which one should be selected and decay
440// the nucleus.
441//
442      G4double totalProb1      = 0.0;
443      G4double probEvapType[6] = {0.0};
444      for (G4int ich=0; ich<nChannels; ich++)
445      {
446        totalProb1      += theChannels[ich]->GetEmissionProbability();
447        probEvapType[ich]  = totalProb1 / totalProb;
448      }
449      G4double xi = G4UniformRand();
450      G4int i     = 0;
451      for (i=0; i<nChannels; i++)
452        if (xi < probEvapType[i]) break;
453      if (i > nChannels) i = nChannels - 1;
454      G4FragmentVector *evaporationResult = theChannels[i]->
455        BreakUp(*intermediateNucleus);
456      fragmentVector->push_back((*evaporationResult)[0]);
457      *intermediateNucleus = *(*evaporationResult)[1];
458      delete evaporationResult->back();
459      delete evaporationResult;
460      evapType.erase(iters[i]);
461    }
462    else
463    {
464//
465//
466// Probability for further evaporation is nil so have to escape from this
467// routine and set the energies of the secondaries to 10eV.
468//
469      evaporate = false;
470    }
471  }
472 
473  return;
474}
475////////////////////////////////////////////////////////////////////////////////
476//
477void G4WilsonAblationModel::SelectSecondariesByDefault (G4ThreeVector boost)
478{
479  for (unsigned i=0; i<evapType.size(); i++)
480  {
481    G4ParticleDefinition *type = fragType[i];
482    G4double mass              = type->GetPDGMass();
483    G4double e                 = mass + 10.0*eV;
484    G4double p                 = std::sqrt(e*e-mass*mass);
485    G4double costheta          = 2.0*G4UniformRand() - 1.0;
486    G4double sintheta          = std::sqrt((1.0 - costheta)*(1.0 + costheta));
487    G4double phi               = twopi * G4UniformRand() * rad;
488    G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
489    G4LorentzVector lorentzVector = G4LorentzVector(direction*p, e);
490    lorentzVector.boost(-boost);
491    G4Fragment *fragment          = 
492      new G4Fragment(lorentzVector, type);
493    fragmentVector->push_back(fragment);
494  }
495}
496////////////////////////////////////////////////////////////////////////////////
497//
498void G4WilsonAblationModel::PrintWelcomeMessage ()
499{
500  G4cout <<G4endl;
501  G4cout <<" *****************************************************************"
502         <<G4endl;
503  G4cout <<" Nuclear ablation model for nuclear-nuclear interactions activated"
504         <<G4endl;
505  G4cout <<" (Written by QinetiQ Ltd for the European Space Agency)"
506         <<G4endl;
507  G4cout <<" *****************************************************************"
508         <<G4endl;
509  G4cout << G4endl;
510
511  return;
512}
513////////////////////////////////////////////////////////////////////////////////
514//
Note: See TracBrowser for help on using the repository browser.