| 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 // Physics model class G4NeutronFissionVI
29 // Created: 03 October 2023
30 // Author V.Ivanchenko
31 //
32
33 #include "G4NeutronFissionVI.hh"
34 #include "G4HadronicInteractionRegistry.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "G4ParticleDefinition.hh"
37 #include "G4Fragment.hh"
38 #include "G4FragmentVector.hh"
39 #include "G4NucleiProperties.hh"
40 #include "G4VEvaporation.hh"
41 #include "G4VEvaporationChannel.hh"
42 #include "G4VPreCompoundModel.hh"
43 #include "G4ExcitationHandler.hh"
44 #include "G4PhotonEvaporation.hh"
45 #include "G4DynamicParticle.hh"
46 #include "G4ParticleHPManager.hh"
47 #include "G4ParticleTable.hh"
48 #include "G4IonTable.hh"
49 #include "G4Electron.hh"
50 #include "G4Deuteron.hh"
51 #include "G4Triton.hh"
52 #include "G4He3.hh"
53 #include "G4Alpha.hh"
54 #include "Randomize.hh"
55 #include "G4RandomDirection.hh"
56 #include "G4HadronicParameters.hh"
57 #include "G4PhysicsModelCatalog.hh"
58
59 G4NeutronFissionVI::G4NeutronFissionVI()
60 : G4HadronicInteraction("nFissionVI"),
61 fManagerHP(G4ParticleHPManager::GetInstance()),
62 minExcitation(0.1*CLHEP::keV),
63 emaxT(fManagerHP->GetMaxEnergyDoppler()),
64 lab4mom(0.,0.,0.,0.)
65 {
66 SetMinEnergy( 0.0*CLHEP::GeV );
67 SetMaxEnergy( G4HadronicParameters::Instance()->GetMaxEnergy() );
68 theTableOfIons = G4ParticleTable::GetParticleTable()->GetIonTable();
69 }
70
71 G4NeutronFissionVI::~G4NeutronFissionVI()
72 {
73 if (fLocalHandler) { delete fHandler; }
74 }
75
76 void G4NeutronFissionVI::InitialiseModel()
77 {
78 if (fFission != nullptr && fHandler != nullptr) { return; }
79 G4HadronicInteraction* p =
80 G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
81 if (nullptr != p) {
82 fHandler = (static_cast<G4VPreCompoundModel*>(p))->GetExcitationHandler();
83 }
84 if (nullptr == fHandler) {
85 fHandler = new G4ExcitationHandler();
86 fLocalHandler = true;
87 }
88 fHandler->Initialise();
89 fFission = fHandler->GetEvaporation()->GetFissionChannel();
90
91 G4DeexPrecoParameters* param =
92 G4NuclearLevelData::GetInstance()->GetParameters();
93 minExcitation = param->GetMinExcitation();
94 secID = G4PhysicsModelCatalog::GetModelID("model_" + GetModelName());
95 }
96
97 G4HadFinalState* G4NeutronFissionVI::ApplyYourself(
98 const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
99 {
100 theParticleChange.Clear();
101 theParticleChange.SetStatusChange(stopAndKill);
102 G4double T = aTrack.GetMaterial()->GetTemperature();
103
104 G4int A = theNucleus.GetA_asInt();
105 G4int Z = theNucleus.GetZ_asInt();
106
107 G4double time = aTrack.GetGlobalTime();
108 G4double ekin = aTrack.GetKineticEnergy();
109
110 // Create initial state
111 G4double mass = G4NucleiProperties::GetNuclearMass(A, Z);
112
113 // no Doppler broading
114 G4double factT = T/CLHEP::STP_Temperature;
115
116 if (ekin >= emaxT*factT || fManagerHP->GetNeglectDoppler()) {
117 lab4mom.set(0.,0.,0.,mass);
118
119 } else {
120 G4double lambda = 1.0/(CLHEP::k_Boltzmann*T);
121 G4double erand = G4RandGamma::shoot(2.0, lambda);
122 auto mom = G4RandomDirection()*std::sqrt(2*mass*erand);
123 lab4mom.set(mom.x(), mom.y(), mom.z(), mass + erand);
124 }
125
126 lab4mom += aTrack.Get4Momentum();
127
128 G4double M = lab4mom.mag();
129 ++A;
130 mass = G4NucleiProperties::GetNuclearMass(A, Z);
131 //G4cout << "Fission start: Z= " << Z << " A= " << A
132 // << " LabM= " << M << " Mcompound= " << mass << G4endl;
133
134 // protection against wrong kinematic
135 if (M < mass) {
136 G4double etot = std::max(mass, lab4mom.e());
137 G4double ptot = std::sqrt((etot - mass)*(etot + mass));
138 G4ThreeVector v = lab4mom.vect().unit();
139 lab4mom.set(v.x()*ptot,v.y()*ptot,v.z()*ptot,etot);
140 }
141
142 G4Fragment* aFragment = new G4Fragment(A, Z, lab4mom);
143
144 if (verboseLevel > 1) {
145 G4cout << "G4NeutronFissionVI::ApplyYourself initial G4Fragmet:"
146 << G4endl;
147 G4cout << aFragment << G4endl;
148 }
149
150 //
151 // Sample final state
152 //
153 fFission->GetEmissionProbability(aFragment);
154 G4Fragment* frag = fFission->EmittedFragment(aFragment);
155 G4ReactionProductVector* final = fHandler->BreakItUp(*aFragment);
156 if (nullptr != frag) {
157 G4ReactionProductVector* v = fHandler->BreakItUp(*frag);
158 for (auto & p : *v) {
159 final->push_back(p);
160 }
161 delete v;
162 delete frag;
163 }
164
165 if (verboseLevel > 1) {
166 G4cout << "G4NeutronFissionVI: " << final->size()
167 << " final particle secID= " << secID << G4endl;
168 }
169 for (auto const & ptr : *final) {
170 G4double etot = ptr->GetTotalEnergy();
171 const G4ParticleDefinition* theDef = ptr->GetDefinition();
172 ekin = std::max(0.0, etot - theDef->GetPDGMass());
173 if (verboseLevel > 1) {
174 G4cout << theDef->GetParticleName()
175 << " Ekin(MeV)= " << ekin/MeV
176 << " p: " << ptr->GetMomentum()
177 << G4endl;
178 }
179 G4HadSecondary* news = new G4HadSecondary(
180 new G4DynamicParticle(theDef, ptr->GetMomentum().unit(), ekin));
181 G4double timeF = std::max(ptr->GetFormationTime(), 0.0);
182 news->SetTime(time + timeF);
183 news->SetCreatorModelID(secID);
184 theParticleChange.AddSecondary(*news);
185 delete news;
186 delete ptr;
187 }
188 delete final;
189 delete aFragment;
190
191 //G4cout << "Fission done" << G4endl;
192 return &theParticleChange;
193 }
194
195