Geant4_10
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 76888 2013-11-18 13:00:28Z 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 "G4PhotonEvaporation.hh"
44 #include "G4DynamicParticle.hh"
45 #include "G4ParticleDefinition.hh"
46 #include "G4ParticleTable.hh"
47 #include "G4IonTable.hh"
48 #include "G4Deuteron.hh"
49 #include "G4Triton.hh"
50 #include "G4He3.hh"
51 #include "G4Alpha.hh"
52 
54  : G4HadronicInteraction("nRadCapture")
55 {
56  lowestEnergyLimit = 0.1*eV;
57  minExcitation = 1*keV;
58  SetMinEnergy( 0.0*GeV );
59  SetMaxEnergy( 100.*TeV );
60  photonEvaporation = new G4PhotonEvaporation();
61  //photonEvaporation = 0;
62  theTableOfIons = G4ParticleTable::GetParticleTable()->GetIonTable();
63 }
64 
66 {
67  delete photonEvaporation;
68 }
69 
71  const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
72 {
75 
76  G4int A = theNucleus.GetA_asInt();
77  G4int Z = theNucleus.GetZ_asInt();
78 
79  G4double time = aTrack.GetGlobalTime();
80 
81  // Create initial state
83  G4LorentzVector lv0(0.0,0.0,0.0,m1);
84  G4LorentzVector lv1 = aTrack.Get4Momentum() + lv0;
85 
86  // simplified method of 1 gamma emission
87  if(A <= 3) {
88 
89  G4ThreeVector bst = lv1.boostVector();
90  G4double M = lv1.mag();
91 
92  ++A;
94  if(M - mass <= lowestEnergyLimit) {
95  return &theParticleChange;
96  }
97 
98  if (verboseLevel > 1) {
99  G4cout << "G4NeutronRadCapture::DoIt: Eini(MeV)="
100  << aTrack.GetKineticEnergy()/MeV << " Eexc(MeV)= "
101  << (M - mass)/MeV
102  << " Z= " << Z << " A= " << A << G4endl;
103  }
104  G4double e1 = (M - mass)*(M + mass)/(2*M);
105  G4double cost = 2.0*G4UniformRand() - 1.0;
106  if(cost > 1.0) {cost = 1.0;}
107  else if(cost < -1.0) {cost = -1.0;}
108  G4double sint = std::sqrt((1. - cost)*(1.0 + cost));
109  G4double phi = G4UniformRand()*CLHEP::twopi;
110  G4LorentzVector lv2(e1*sint*std::cos(phi),e1*sint*std::sin(phi),e1*cost,e1);
111  lv2.boost(bst);
112  G4HadSecondary* news =
114  news->SetTime(time);
116  delete news;
117 
118  G4ParticleDefinition* theDef = 0;
119 
120  lv1 -= lv2;
121  if (Z == 1 && A == 2) {theDef = G4Deuteron::Deuteron();}
122  else if (Z == 1 && A == 3) {theDef = G4Triton::Triton();}
123  else if (Z == 2 && A == 3) {theDef = G4He3::He3();}
124  else if (Z == 2 && A == 4) {theDef = G4Alpha::Alpha();}
125  else { theDef = theTableOfIons->GetIon(Z,A,0); }
126 
127  if (verboseLevel > 1) {
128  G4cout << "Gamma 4-mom: " << lv2 << " "
129  << theDef->GetParticleName() << " " << lv1 << G4endl;
130  }
131  if(theDef) {
132  news = new G4HadSecondary(new G4DynamicParticle(theDef, lv1));
133  news->SetTime(time);
135  delete news;
136  }
137 
138  // Use photon evaporation
139  } else {
140 
141  G4Fragment* aFragment = new G4Fragment(A+1, Z, lv1);
142 
143  if (verboseLevel > 1) {
144  G4cout << "G4NeutronRadCapture::ApplyYourself initial G4Fragmet:" << G4endl;
145  G4cout << aFragment << G4endl;
146  }
147 
148  //
149  // Sample final state
150  //
151  G4FragmentVector* fv = photonEvaporation->BreakUpFragment(aFragment);
152  if(!fv) { fv = new G4FragmentVector(); }
153  fv->push_back(aFragment);
154  size_t n = fv->size();
155 
156  if (verboseLevel > 1) {
157  G4cout << "G4NeutronRadCapture: " << n << " final particle" << G4endl;
158  }
159  for(size_t i=0; i<n; ++i) {
160 
161  G4Fragment* f = (*fv)[i];
162  G4double etot = f->GetMomentum().e();
163 
164  Z = f->GetZ_asInt();
165  A = f->GetA_asInt();
166 
167  G4ParticleDefinition* theDef = 0;
168  if(0 == Z && 0 == A) {theDef = f->GetParticleDefinition();}
169  else if (Z == 1 && A == 2) {theDef = G4Deuteron::Deuteron();}
170  else if (Z == 1 && A == 3) {theDef = G4Triton::Triton();}
171  else if (Z == 2 && A == 3) {theDef = G4He3::He3();}
172  else if (Z == 2 && A == 4) {theDef = G4Alpha::Alpha();}
173  else {
174  G4double eexc = f->GetExcitationEnergy();
175  G4double excitation = eexc;
176  G4int level = 0;
177  theDef = theTableOfIons->GetIon(Z, A, level);
178  /*
179  G4cout << "### Find ion Z= " << theFragmentZ << " A= " << theFragmentA
180  << " Eexc(MeV)= " << excitation/MeV << " "
181  << theKindOfFragment << G4endl;
182  */
183  // production of an isomer
184  if(eexc > minExcitation) {
185  G4double elevel1 = 0.0;
186  G4double elevel2 = 0.0;
188  for(level=1; level<9; ++level) {
189  ion = theTableOfIons->GetIon(Z, A, level);
190  //G4cout << level << " " << ion << G4endl;
191  if(ion) {
192  G4Ions* ip = dynamic_cast<G4Ions*>(ion);
193  if(ip) {
194  elevel2 = ip->GetExcitationEnergy();
195  //G4cout<<" Level "<<level<<" E(MeV)= "<<elevel2/MeV<<G4endl;
196  // close level
197  if(std::fabs(eexc - elevel2) < minExcitation) {
198  excitation = eexc - elevel2;
199  theDef = ion;
200  break;
201  // previous level was closer
202  } else if(elevel2 - eexc >= eexc - elevel1) {
203  excitation = eexc - elevel1;
204  break;
205  // will check next level and save current
206  } else {
207  theDef = ion;
208  excitation = eexc - elevel2;
209  elevel1 = elevel2;
210  }
211  }
212  } else {
213  break;
214  }
215  }
216  }
217  // correction of total energy for ground state isotopes
218  etot += excitation;
219  etot -= theDef->GetPDGMass();
220  if(etot < 0.0) { etot = 0.0; }
221  }
222  if (verboseLevel > 1) {
223  G4cout << i << ". " << theDef->GetParticleName()
224  << " Ekin(MeV)= " << etot/MeV
225  << " p: " << f->GetMomentum().vect()
226  << G4endl;
227  }
228  if(theDef) {
229  G4HadSecondary* news =
230  new G4HadSecondary(new G4DynamicParticle(theDef,
231  f->GetMomentum().vect().unit(),
232  etot));
233  G4double timeF = f->GetCreationTime();
234  if(timeF < 0.0) { timeF = 0.0; }
235  news->SetTime(time + timeF);
237  delete news;
238  }
239  delete f;
240  }
241  delete fv;
242  }
243  return &theParticleChange;
244 }
245 
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
Hep3Vector boostVector() const
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double GetExcitationEnergy() const
Definition: G4Ions.hh:113
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:449
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
TFile f
Definition: plotHisto.C:6
void SetStatusChange(G4HadFinalStateStatus aS)
void SetMinEnergy(G4double anEnergy)
Char_t n[5]
G4IonTable * GetIonTable() const
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
G4int GetA_asInt() const
Definition: G4Fragment.hh:218
Definition: G4Ions.hh:51
G4double GetCreationTime() const
Definition: G4Fragment.hh:378
Float_t Z
Definition: plot.C:39
double mag() const
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:251
HepLorentzVector & boost(double, double, double)
G4double GetKineticEnergy() const
G4double GetGlobalTime() const
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:64
static G4Triton * Triton()
Definition: G4Triton.cc:95
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
const G4LorentzVector & Get4Momentum() const
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
void SetTime(G4double aT)
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
Hep3Vector unit() const
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4int GetZ_asInt() const
Definition: G4Fragment.hh:223
G4ParticleDefinition * GetParticleDefinition() const
Definition: G4Fragment.hh:368
void SetMaxEnergy(const G4double anEnergy)
#define G4endl
Definition: G4ios.hh:61
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
virtual G4FragmentVector * BreakUpFragment(G4Fragment *theNucleus)
double G4double
Definition: G4Types.hh:76
static G4He3 * He3()
Definition: G4He3.cc:94
void AddSecondary(G4DynamicParticle *aP)
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:235