| Geant4 Cross Reference |
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 //
27 // -------------------------------------------------------------------
28 // GEANT 4 class file
29 //
30 // CERN, Geneva, Switzerland
31 //
32 // File name: G4UnstableFragmentBreakUp
33 //
34 // Author: V.Ivanchenko
35 //
36 // Creation date: 11 May 2010
37 //
38 //Modifications:
39 //
40 // -------------------------------------------------------------------
41 //
42
43 #include "G4UnstableFragmentBreakUp.hh"
44 #include "Randomize.hh"
45 #include "G4RandomDirection.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4LorentzVector.hh"
48 #include "G4Fragment.hh"
49 #include "G4FragmentVector.hh"
50 #include "G4NucleiProperties.hh"
51 #include "G4NuclearLevelData.hh"
52 #include "G4LevelManager.hh"
53 #include "Randomize.hh"
54 #include "G4PhysicsModelCatalog.hh"
55
56 const G4int G4UnstableFragmentBreakUp::Zfr[] = {0, 1, 1, 1, 2, 2};
57 const G4int G4UnstableFragmentBreakUp::Afr[] = {1, 1, 2, 3, 3, 4};
58
59 G4UnstableFragmentBreakUp::G4UnstableFragmentBreakUp() : fVerbose(1), fSecID(-1)
60 {
61 fLevelData = G4NuclearLevelData::GetInstance();
62 for(G4int i=0; i<6; ++i) {
63 masses[i] = G4NucleiProperties::GetNuclearMass(Afr[i], Zfr[i]);
64 }
65 fSecID = G4PhysicsModelCatalog::GetModelID("model_G4UnstableFragmentBreakUp");
66 }
67
68 G4UnstableFragmentBreakUp::~G4UnstableFragmentBreakUp()
69 {}
70
71 G4bool G4UnstableFragmentBreakUp::BreakUpChain(G4FragmentVector* results,
72 G4Fragment* nucleus)
73 {
74 //G4cout << "G4UnstableFragmentBreakUp::EmittedFragment" << G4endl;
75 G4Fragment* frag = nullptr;
76
77 G4int Z = nucleus->GetZ_asInt();
78 G4int A = nucleus->GetA_asInt();
79
80 G4LorentzVector lv = nucleus->GetMomentum();
81 G4double time = nucleus->GetCreationTime();
82
83 G4double mass1(0.0), mass2(0.0);
84
85 // look for the decay channel with normal masses
86 // without Coulomb barrier and paring corrections
87 // 1 - recoil, 2 - emitted light ion
88 if(fVerbose > 1) {
89 G4cout << "#Unstable decay " << " Z= " << Z << " A= " << A
90 << " Eex(MeV)= " << nucleus->GetExcitationEnergy() << G4endl;
91 }
92 const G4double tolerance = 10*CLHEP::eV;
93 const G4double dmlimit = 0.005*CLHEP::MeV;
94 G4double mass = lv.mag();
95 G4double exca = -1000.0;
96 G4bool isChannel = false;
97 G4int idx = -1;
98 for(G4int i=0; i<6; ++i) {
99 G4int Zres = Z - Zfr[i];
100 G4int Ares = A - Afr[i];
101 if(Zres >= 0 && Ares >= Zres && Ares >= Afr[i]) {
102 if(Ares <= 4) {
103 for(G4int j=0; j<6; ++j) {
104 if(Zres == Zfr[j] && Ares == Afr[j]) {
105 G4double delm = mass - masses[i] - masses[j];
106 /*
107 G4cout << "i=" << i << " j=" << j << " Zres=" << Zres
108 << " Ares=" << Ares << " delm=" << delm << G4endl;
109 */
110 if(delm > exca) {
111 mass2 = masses[i]; // emitted
112 mass1 = masses[j]; // recoil
113 exca = delm;
114 idx = i;
115 if(delm > 0.0) {
116 isChannel = true;
117 break;
118 }
119 }
120 }
121 }
122 }
123 if(isChannel) { break; }
124 // no simple channel
125 G4double mres = G4NucleiProperties::GetNuclearMass(Ares, Zres);
126 G4double e = mass - mres - masses[i];
127 // select excited state
128 //const G4LevelManager* lman = fLevelData->GetLevelManager(Zres, Ares);
129 /*
130 G4cout << "i=" << i << " Zres=" << Zres
131 << " Ares=" << Ares << " delm=" << e << G4endl;
132 */
133 if(e >= exca) {
134 mass2 = masses[i];
135 mass1 = (Ares > 4 && e > 0.0) ? mres + e*G4UniformRand() : mres;
136 exca = e;
137 idx = i;
138 if(e > 0.0) {
139 isChannel = true;
140 break;
141 }
142 }
143 }
144 }
145
146 G4double massmin = mass1 + mass2;
147 if(fVerbose > 1) {
148 G4cout << "isChannel:" << isChannel << " idx=" << idx << " Zfr=" << Zfr[idx]
149 << " Arf=" << Afr[idx] << " delm=" << mass - massmin << G4endl;
150 }
151 if(!isChannel || mass < massmin) {
152 if(mass + dmlimit < massmin) { return false; }
153 if(fVerbose > 1) {
154 G4cout << "#Unstable decay correction: Z= " << Z << " A= " << A
155 << " idx= " << idx
156 << " deltaM(MeV)= " << mass - massmin
157 << G4endl;
158 }
159 mass = massmin;
160 G4double e = std::max(lv.e(), mass + tolerance);
161 G4double mom = std::sqrt((e - mass)*(e + mass));
162 G4ThreeVector dir = lv.vect().unit();
163 lv.set(dir*mom, e);
164 }
165
166 // compute energy of light fragment
167 G4double e2 = 0.5*((mass - mass1)*(mass + mass1) + mass2*mass2)/mass;
168 e2 = std::max(e2, mass2);
169 G4double mom = std::sqrt((e2 - mass2)*(e2 + mass2));
170
171 // sample decay
172 G4ThreeVector bst = lv.boostVector();
173 G4ThreeVector v = G4RandomDirection();
174 G4LorentzVector mom2 = G4LorentzVector(v*mom, e2);
175 mom2.boost(bst);
176 frag = new G4Fragment(Afr[idx], Zfr[idx], mom2);
177 frag->SetCreationTime(time);
178 frag->SetCreatorModelID(fSecID);
179 results->push_back(frag);
180
181 // residual
182 lv -= mom2;
183 Z -= Zfr[idx];
184 A -= Afr[idx];
185
186 nucleus->SetZandA_asInt(Z, A);
187 nucleus->SetMomentum(lv);
188 nucleus->SetCreatorModelID(fSecID);
189 return true;
190 }
191
192 G4double G4UnstableFragmentBreakUp::GetEmissionProbability(G4Fragment*)
193 {
194 return 0.0;
195 }
196