Geant4  10.03
G4UnstableFragmentBreakUp.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: G4UnstableFragmentBreakUp.cc 99812 2016-10-06 08:49:27Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 // GEANT 4 class file
30 //
31 // CERN, Geneva, Switzerland
32 //
33 // File name: G4UnstableFragmentBreakUp
34 //
35 // Author: V.Ivanchenko
36 //
37 // Creation date: 11 May 2010
38 //
39 //Modifications:
40 //
41 // -------------------------------------------------------------------
42 //
43 
45 #include "Randomize.hh"
46 #include "G4RandomDirection.hh"
47 #include "G4PhysicalConstants.hh"
48 #include "G4LorentzVector.hh"
49 #include "G4Fragment.hh"
50 #include "G4FragmentVector.hh"
51 #include "G4NucleiProperties.hh"
52 #include "G4NuclearLevelData.hh"
53 #include "G4LevelManager.hh"
54 
55 const G4int G4UnstableFragmentBreakUp::Zfr[] = {0, 1, 1, 1, 2, 2};
56 const G4int G4UnstableFragmentBreakUp::Afr[] = {1, 1, 2, 3, 3, 4};
57 
59 {
61  for(G4int i=0; i<6; ++i) {
63  }
64 }
65 
67 {}
68 
70  G4Fragment* nucleus)
71 {
72  //G4cout << "G4UnstableFragmentBreakUp::EmittedFragment" << G4endl;
73  G4Fragment* frag = nullptr;
74 
75  G4int Z = nucleus->GetZ_asInt();
76  G4int A = nucleus->GetA_asInt();
77 
78  // if the isotope is in the database it is not exotic
79  // so, cannot be handled by this class
80  if(fLevelData->GetLevelManager(Z, A)) { return false; }
81 
82  G4LorentzVector lv = nucleus->GetMomentum();
83  G4double time = nucleus->GetCreationTime();
84 
85  G4double mass1(0.0), mass2(0.0);
86 
87  G4double mass = lv.mag();
88  G4bool done = false;
89 
90  G4int Amax = A;
91  G4int Zres = Z;
92  G4int Ares = A;
93 
94  for(G4int k=0; k<Amax; ++k) {
95  // look for the decay channel with normal masses
96  // without Coulomb barrier and paring corrections
97 
98  // 1 - recoil, 2 - emitted light ion
99 
100  //G4cout << "Unstable decay #" << k << " Z= " << Z << " A= " << A << G4endl;
101 
102  G4double ekin = 0.0;
103  G4int index = -1;
104  for(G4int i=0; i<6; ++i) {
105  Zres = Z - Zfr[i];
106  Ares = A - Afr[i];
107  for(G4int j=0; j<6; ++j) {
108  if(Zres == Zfr[j] && Ares == Afr[j]) {
109  /*
110  G4cout << "i= " << i << " j= " << j << " Zres= " << Zres
111  << " Ares= " << Ares << " dm= " << mass - masses[i] - masses[j]
112  << G4endl;
113  */
114  if(mass >= masses[i] + masses[j]) {
115  mass2 = masses[i];
116  mass1 = masses[j];
117  ekin = 0.0;
118  done = true;
119  break;
120  }
121  }
122  }
123  if(done) {
124  index = i;
125  break;
126  }
127  if(Zres >= 0 && Ares >= Zres && Ares > 0) {
129  G4double e = mass - m - masses[i];
130  if(e > ekin) {
131  mass2 = masses[i];
132  mass1 = m;
133  ekin = e;
134  index = i;
135  if(fLevelData->GetLevelManager(Z, A)) {
136  done = true;
137  break;
138  }
139  }
140  }
141  }
142 
143  // no decay channel - assume that primary mass is biased
144  // only energy will be conserved
145  if(index < 0) {
146  for(G4int i=0; i<6; ++i) {
147  Zres = Z - Zfr[i];
148  Ares = A - Afr[i];
149  if(Zres >= 0 && Ares >= Zres && Ares > 0) {
150 
151  G4double m = proton_mass_c2*Zres + neutron_mass_c2*(Ares - Zres);
152  G4double e = mass - m - masses[i];
153  if(e > ekin) {
154  mass2 = masses[i];
155  mass1 = m;
156  ekin = e;
157  index = i;
158  }
159  }
160  }
161  }
162  if(index < 0) {
163 
164  // further decay impossible
165  // if no one decay sampled do not update primary
166  if(0 == k) { return false; }
167  else { break; }
168  }
169 
170  // useful to left max excitation for the residual
171  mass1 += ekin;
172 
173  // compute energy of light fragment
174  G4double e2 = 0.5*((mass - mass1)*(mass + mass1) + mass2*mass2)/mass;
175  e2 = std::max(e2, mass2);
176  G4double mom = std::sqrt((e2 - mass2)*(e2 + mass2));
177 
178  // sample decay
179  G4ThreeVector bst = lv.boostVector();
181  G4LorentzVector mom2 = G4LorentzVector(v*mom, e2);
182  mom2.boost(bst);
183  frag = new G4Fragment(Afr[index], Zfr[index], mom2);
184  frag->SetCreationTime(time);
185  results->push_back(frag);
186 
187  // residual
188  lv -= mom2;
189  Z -= Zfr[index];
190  A -= Afr[index];
191 
192  mass = lv.mag();
193 
194  if(done) { break; }
195  }
196 
197  // updated primary
198  nucleus->SetZandA_asInt(Z, A);
199  nucleus->SetMomentum(lv);
200  return true;
201 }
202 
204 {
205  return 0.0;
206 }
static G4double GetNuclearMass(const G4double A, const G4double Z)
CLHEP::Hep3Vector G4ThreeVector
const G4LevelManager * GetLevelManager(G4int Z, G4int A)
G4ThreeVector G4RandomDirection()
virtual G4bool BreakUpChain(G4FragmentVector *, G4Fragment *) final
int G4int
Definition: G4Types.hh:78
double A(double temperature)
G4int GetA_asInt() const
Definition: G4Fragment.hh:256
static constexpr double m
Definition: G4SIunits.hh:129
G4double GetCreationTime() const
Definition: G4Fragment.hh:420
bool G4bool
Definition: G4Types.hh:79
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:289
void SetMomentum(const G4LorentzVector &value)
Definition: G4Fragment.hh:294
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:63
virtual G4double GetEmissionProbability(G4Fragment *fragment) final
void SetCreationTime(G4double time)
Definition: G4Fragment.hh:425
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void SetZandA_asInt(G4int Znew, G4int Anew)
Definition: G4Fragment.hh:266
G4int GetZ_asInt() const
Definition: G4Fragment.hh:261
double G4double
Definition: G4Types.hh:76
static G4NuclearLevelData * GetInstance()
CLHEP::HepLorentzVector G4LorentzVector