Geant4  10.02.p01
G4NeutronRadCapture.cc
Go to the documentation of this file.
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: G4NeutronRadCapture.cc 88191 2015-02-02 17:27:37Z gcosmo $
27 //
28 //
29 // Physics model class G4NeutronRadCapture
30 // Created: 31 August 2009
31 // Author V.Ivanchenko
32 //
33 // Modified:
34 // 09.09.2010 V.Ivanchenko added usage of G4PhotonEvaporation
35 //
36 
37 #include "G4NeutronRadCapture.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4ParticleDefinition.hh"
40 #include "G4Fragment.hh"
41 #include "G4FragmentVector.hh"
42 #include "G4NucleiProperties.hh"
43 #include "G4VEvaporationChannel.hh"
44 #include "G4PhotonEvaporation.hh"
46 #include "G4DynamicParticle.hh"
47 #include "G4ParticleTable.hh"
48 #include "G4IonTable.hh"
49 #include "G4Deuteron.hh"
50 #include "G4Triton.hh"
51 #include "G4He3.hh"
52 #include "G4Alpha.hh"
53 
55  : G4HadronicInteraction("nRadCapture"),
56  lab4mom(0.,0.,0.,0.)
57 {
58  lowestEnergyLimit = 10*eV;
59  minExcitation = 1*keV;
60  SetMinEnergy( 0.0*GeV );
61  SetMaxEnergy( 100.*TeV );
62 
63  char* env = getenv("G4UsePhotonEvaporationOLD");
64  if(!env) { photonEvaporation = new G4PhotonEvaporation(); }
67 
69 }
70 
72 {
73  delete photonEvaporation;
74 }
75 
77  const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
78 {
81 
82  G4int A = theNucleus.GetA_asInt();
83  G4int Z = theNucleus.GetZ_asInt();
84 
85  G4double time = aTrack.GetGlobalTime();
86 
87  // Create initial state
88  lab4mom.set(0.,0.,0.,G4NucleiProperties::GetNuclearMass(A, Z));
89  lab4mom += aTrack.Get4Momentum();
90 
91  G4double M = lab4mom.mag();
92  ++A;
94  //G4cout << "Capture start: Z= " << Z << " A= " << A
95  // << " LabM= " << M << " Mcompound= " << mass << G4endl;
96 
97  // simplified method of 1 gamma emission
98  if(A <= 4) {
99 
100  G4ThreeVector bst = lab4mom.boostVector();
101 
102  if(M - mass <= lowestEnergyLimit) {
103  return &theParticleChange;
104  }
105 
106  if (verboseLevel > 1) {
107  G4cout << "G4NeutronRadCapture::DoIt: Eini(MeV)="
108  << aTrack.GetKineticEnergy()/MeV << " Eexc(MeV)= "
109  << (M - mass)/MeV
110  << " Z= " << Z << " A= " << A << G4endl;
111  }
112  G4double e1 = (M - mass)*(M + mass)/(2*M);
113 
114  G4double cost = 2.0*G4UniformRand() - 1.0;
115  if(cost > 1.0) {cost = 1.0;}
116  else if(cost < -1.0) {cost = -1.0;}
117  G4double sint = std::sqrt((1. - cost)*(1.0 + cost));
119 
120  G4LorentzVector lv2(e1*sint*std::cos(phi),e1*sint*std::sin(phi),
121  e1*cost,e1);
122  lv2.boost(bst);
123  G4HadSecondary* news =
125  news->SetTime(time);
127  delete news;
128 
129  const G4ParticleDefinition* theDef = 0;
130 
131  lab4mom -= lv2;
132  if (Z == 1 && A == 2) {theDef = G4Deuteron::Deuteron();}
133  else if (Z == 1 && A == 3) {theDef = G4Triton::Triton();}
134  else if (Z == 2 && A == 3) {theDef = G4He3::He3();}
135  else if (Z == 2 && A == 4) {theDef = G4Alpha::Alpha();}
136  else { theDef = theTableOfIons->GetIon(Z,A,0.0); }
137 
138  if (verboseLevel > 1) {
139  G4cout << "Gamma 4-mom: " << lv2 << " "
140  << theDef->GetParticleName() << " " << lab4mom << G4endl;
141  }
142  if(theDef) {
143  news = new G4HadSecondary(new G4DynamicParticle(theDef, lab4mom));
144  news->SetTime(time);
146  delete news;
147  }
148 
149  // Use photon evaporation
150  } else {
151 
152  // protection against wrong kinematic
153  if(M < mass) {
154  G4double etot = std::max(mass, lab4mom.e());
155  G4double ptot = std::sqrt((etot - mass)*(etot + mass));
156  G4ThreeVector v = lab4mom.vect().unit();
157  lab4mom.set(v.x()*ptot,v.y()*ptot,v.z()*ptot,etot);
158  }
159 
160  G4Fragment* aFragment = new G4Fragment(A, Z, lab4mom);
161 
162  if (verboseLevel > 1) {
163  G4cout << "G4NeutronRadCapture::ApplyYourself initial G4Fragmet:"
164  << G4endl;
165  G4cout << aFragment << G4endl;
166  }
167 
168  //
169  // Sample final state
170  //
172  if(!fv) { fv = new G4FragmentVector(); }
173  fv->push_back(aFragment);
174  size_t n = fv->size();
175 
176  if (verboseLevel > 1) {
177  G4cout << "G4NeutronRadCapture: " << n << " final particle" << G4endl;
178  }
179  for(size_t i=0; i<n; ++i) {
180 
181  G4Fragment* f = (*fv)[i];
182  G4double etot = f->GetMomentum().e();
183 
184  Z = f->GetZ_asInt();
185  A = f->GetA_asInt();
186 
187  const G4ParticleDefinition* theDef;
188  if(0 == Z && 0 == A) {theDef = f->GetParticleDefinition();}
189  else if (Z == 1 && A == 2) {theDef = G4Deuteron::Deuteron();}
190  else if (Z == 1 && A == 3) {theDef = G4Triton::Triton();}
191  else if (Z == 2 && A == 3) {theDef = G4He3::He3();}
192  else if (Z == 2 && A == 4) {theDef = G4Alpha::Alpha();}
193  else {
194  G4double eexc = f->GetExcitationEnergy();
195  if(eexc <= minExcitation) { eexc = 0.0; }
196  theDef = theTableOfIons->GetIon(Z, A, eexc);
197  /*
198  G4cout << "### Find ion Z= " << Z << " A= " << A
199  << " Eexc(MeV)= " << eexc/MeV << " "
200  << theDef << G4endl;
201  */
202  }
203  G4double ekin = std::max(0.0,etot - theDef->GetPDGMass());
204  if (verboseLevel > 1) {
205  G4cout << i << ". " << theDef->GetParticleName()
206  << " Ekin(MeV)= " << etot/MeV
207  << " p: " << f->GetMomentum().vect()
208  << G4endl;
209  }
210  G4HadSecondary* news = new G4HadSecondary(
211  new G4DynamicParticle(theDef,
212  f->GetMomentum().vect().unit(),
213  ekin));
214  G4double timeF = f->GetCreationTime();
215  if(timeF < 0.0) { timeF = 0.0; }
216  news->SetTime(time + timeF);
218  delete news;
219  delete f;
220  }
221  delete fv;
222  }
223  return &theParticleChange;
224 }
225 
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
static const double MeV
Definition: G4SIunits.hh:211
static G4double GetNuclearMass(const G4double A, const G4double Z)
virtual G4FragmentVector * BreakUpFragment(G4Fragment *theNucleus)
CLHEP::Hep3Vector G4ThreeVector
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:491
const G4ParticleDefinition * GetParticleDefinition() const
Definition: G4Fragment.hh:395
virtual void SetICM(G4bool)
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
G4VEvaporationChannel * photonEvaporation
void SetStatusChange(G4HadFinalStateStatus aS)
void SetMinEnergy(G4double anEnergy)
G4IonTable * GetIonTable() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
G4int GetA_asInt() const
Definition: G4Fragment.hh:251
G4double GetCreationTime() const
Definition: G4Fragment.hh:405
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:284
G4double GetKineticEnergy() const
static const double twopi
Definition: G4SIunits.hh:75
G4double GetGlobalTime() const
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:63
static G4Triton * Triton()
Definition: G4Triton.cc:95
static const double GeV
Definition: G4SIunits.hh:214
const G4int n
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
const G4LorentzVector & Get4Momentum() const
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
static const G4double e1
void SetTime(G4double aT)
static const double eV
Definition: G4SIunits.hh:212
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4int GetZ_asInt() const
Definition: G4Fragment.hh:256
void SetMaxEnergy(const G4double anEnergy)
#define G4endl
Definition: G4ios.hh:61
static const double TeV
Definition: G4SIunits.hh:215
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
static const double keV
Definition: G4SIunits.hh:213
double G4double
Definition: G4Types.hh:76
static G4He3 * He3()
Definition: G4He3.cc:94
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:268
CLHEP::HepLorentzVector G4LorentzVector