| 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 // G4TheoFSGenerator
28 //
29 // 20110307 M. Kelsey -- Add call to new theTransport->SetPrimaryProjectile()
30 // to provide access to full initial state (for Bertini)
31 // 20110805 M. Kelsey -- Follow change to G4V3DNucleus::GetNucleons()
32
33 #include "G4DynamicParticle.hh"
34 #include "G4TheoFSGenerator.hh"
35 #include "G4ReactionProductVector.hh"
36 #include "G4ReactionProduct.hh"
37 #include "G4IonTable.hh"
38 #include "G4HadronicParameters.hh"
39 #include "G4CRCoalescence.hh"
40 #include "G4HadronicInteractionRegistry.hh"
41
42 G4TheoFSGenerator::G4TheoFSGenerator(const G4String& name)
43 : G4HadronicInteraction(name)
44 , theTransport(nullptr), theHighEnergyGenerator(nullptr)
45 , theQuasielastic(nullptr)
46 , theCosmicCoalescence(nullptr)
47 , theStringModelID(-1)
48 {
49 theParticleChange = new G4HadFinalState;
50 theStringModelID = G4PhysicsModelCatalog::GetModelID( "model_" + name );
51 }
52
53 G4TheoFSGenerator::~G4TheoFSGenerator()
54 {
55 delete theParticleChange;
56 }
57
58 void G4TheoFSGenerator::ModelDescription(std::ostream& outFile) const
59 {
60 outFile << GetModelName() <<" consists of a " << theHighEnergyGenerator->GetModelName()
61 << " string model and a stage to de-excite the excited nuclear fragment.\n<p>"
62 << "The string model simulates the interaction of\n"
63 << "an incident hadron with a nucleus, forming \n"
64 << "excited strings, decays these strings into hadrons,\n"
65 << "and leaves an excited nucleus. \n"
66 << "<p>The string model:\n";
67 theHighEnergyGenerator->ModelDescription(outFile);
68 outFile <<"\n<p>";
69 theTransport->PropagateModelDescription(outFile);
70 }
71
72 G4HadFinalState * G4TheoFSGenerator::ApplyYourself(const G4HadProjectile & thePrimary, G4Nucleus &theNucleus)
73 {
74 // init particle change
75 theParticleChange->Clear();
76 theParticleChange->SetStatusChange(stopAndKill);
77 G4double timePrimary=thePrimary.GetGlobalTime();
78
79 // Temporarily dummy treatment of heavy (charm and bottom) hadron projectiles at low energies.
80 // Cascade models are currently not applicable for heavy hadrons and string models cannot
81 // handle them properly at low energies - let's say safely below ~100 MeV.
82 // In these cases, we return as final state the initial state unchanged.
83 // For most applications, this is a safe simplification, giving that the nearly all
84 // slowly moving charm and bottom hadrons decay before any hadronic interaction can occur.
85 // Note that we prefer not to use G4HadronicParameters::GetMinEnergyTransitionFTF_Cascade()
86 // (typically ~3 GeV) because FTFP works reasonably well below such a value.
87 const G4double energyThresholdForCharmAndBottomHadrons = 100.0*CLHEP::MeV;
88 if ( thePrimary.GetKineticEnergy() < energyThresholdForCharmAndBottomHadrons &&
89 ( thePrimary.GetDefinition()->GetQuarkContent( 4 ) != 0 || // Has charm constituent quark
90 thePrimary.GetDefinition()->GetAntiQuarkContent( 4 ) != 0 || // Has anti-charm constituent anti-quark
91 thePrimary.GetDefinition()->GetQuarkContent( 5 ) != 0 || // Has bottom constituent quark
92 thePrimary.GetDefinition()->GetAntiQuarkContent( 5 ) != 0 ) ) { // Has anti-bottom constituent anti-quark
93 theParticleChange->SetStatusChange( isAlive );
94 theParticleChange->SetEnergyChange( thePrimary.GetKineticEnergy() );
95 theParticleChange->SetMomentumChange( thePrimary.Get4Momentum().vect().unit() );
96 return theParticleChange;
97 }
98
99 // Temporarily dummy treatment of light hypernuclei projectiles at low energies.
100 // Bertini and Binary intra-nuclear cascade models are currently not applicable
101 // for hypernuclei and string models cannot handle them properly at low energies -
102 // let's say safely below ~100 MeV.
103 // In these cases, we return as final state the initial state unchanged.
104 // For most applications, this is a safe simplification, giving that the nearly all
105 // slowly moving hypernuclei decay before any hadronic interaction can occur.
106 // Note that we prefer not to use G4HadronicParameters::GetMinEnergyTransitionFTF_Cascade()
107 // (typically ~3 GeV) because FTFP works reasonably well below such a value.
108 // Note that for anti-hypernuclei, FTF model works fine down to zero kinetic energy,
109 // so there is no need of any extra, dummy treatment.
110 const G4double energyThresholdForHyperNuclei = 100.0*CLHEP::MeV;
111 if ( thePrimary.GetKineticEnergy() < energyThresholdForHyperNuclei &&
112 thePrimary.GetDefinition()->IsHypernucleus() ) {
113 theParticleChange->SetStatusChange( isAlive );
114 theParticleChange->SetEnergyChange( thePrimary.GetKineticEnergy() );
115 theParticleChange->SetMomentumChange( thePrimary.Get4Momentum().vect().unit() );
116 return theParticleChange;
117 }
118
119 // check if models have been registered, and use default, in case this is not true @@
120
121 const G4DynamicParticle aPart(thePrimary.GetDefinition(),thePrimary.Get4Momentum().vect());
122
123 if ( theQuasielastic )
124 {
125 if ( theQuasielastic->GetFraction(theNucleus, aPart) > G4UniformRand() )
126 {
127 //G4cout<<"___G4TheoFSGenerator: before Scatter (1) QE=" << theQuasielastic<<G4endl;
128 G4KineticTrackVector *result= theQuasielastic->Scatter(theNucleus, aPart);
129 //G4cout << "^^G4TheoFSGenerator: after Scatter (1) " << G4endl;
130 if (result)
131 {
132 for(auto & ptr : *result)
133 {
134 G4DynamicParticle * aNew =
135 new G4DynamicParticle(ptr->GetDefinition(),
136 ptr->Get4Momentum().e(),
137 ptr->Get4Momentum().vect());
138 theParticleChange->AddSecondary(aNew, ptr->GetCreatorModelID());
139 delete ptr;
140 }
141 delete result;
142 }
143 else
144 {
145 theParticleChange->SetStatusChange(isAlive);
146 theParticleChange->SetEnergyChange(thePrimary.GetKineticEnergy());
147 theParticleChange->SetMomentumChange(thePrimary.Get4Momentum().vect().unit());
148 }
149 return theParticleChange;
150 }
151 }
152
153 // get result from high energy model
154 G4KineticTrackVector * theInitialResult =
155 theHighEnergyGenerator->Scatter(theNucleus, aPart);
156
157 // Assign the creator model ID
158 for ( auto & ptr : *theInitialResult ) {
159 ptr->SetCreatorModelID( theStringModelID );
160 }
161
162 //#define DEBUG_initial_result
163 #ifdef DEBUG_initial_result
164 G4double E_out(0);
165 G4IonTable * ionTable=G4ParticleTable::GetParticleTable()->GetIonTable();
166 for(auto & ptr : *theInitialResult)
167 {
168 //G4cout << "TheoFS secondary, mom " << ptr->GetDefinition()->GetParticleName()
169 // << " " << ptr->Get4Momentum() << G4endl;
170 E_out += ptr->Get4Momentum().e();
171 }
172 G4double init_mass= ionTable->GetIonMass(theNucleus.GetZ_asInt(),theNucleus.GetA_asInt());
173 G4double init_E=aPart.Get4Momentum().e();
174 // residual nucleus
175
176 const std::vector<G4Nucleon> & thy = theHighEnergyGenerator->GetWoundedNucleus()->GetNucleons();
177
178 G4int resZ(0),resA(0);
179 G4double delta_m(0);
180 for(auto & nuc : thy)
181 {
182 if(nuc.AreYouHit()) {
183 ++resA;
184 if ( nuc.GetDefinition() == G4Proton::Proton() ) {
185 ++resZ;
186 delta_m += CLHEP::proton_mass_c2;
187 } else {
188 delta_m += CLHEP::neutron_mass_c2;
189 }
190 }
191 }
192
193 G4double final_mass(0);
194 if ( theNucleus.GetA_asInt() ) {
195 final_mass=ionTable->GetIonMass(theNucleus.GetZ_asInt()-resZ,theNucleus.GetA_asInt()- resA);
196 }
197 G4double E_excit=init_mass + init_E - final_mass - E_out;
198 G4cout << " Corrected delta mass " << init_mass - final_mass - delta_m << G4endl;
199 G4cout << "initial E, mass = " << init_E << ", " << init_mass << G4endl;
200 G4cout << " final E, mass = " << E_out <<", " << final_mass << " excitation_E " << E_excit << G4endl;
201 #endif
202
203 G4ReactionProductVector * theTransportResult = nullptr;
204
205 G4V3DNucleus* theProjectileNucleus = theHighEnergyGenerator->GetProjectileNucleus();
206 if(theProjectileNucleus == nullptr)
207 {
208 G4int hitCount = 0;
209 const std::vector<G4Nucleon>& they = theHighEnergyGenerator->GetWoundedNucleus()->GetNucleons();
210 for(auto & nuc : they)
211 {
212 if(nuc.AreYouHit()) ++hitCount;
213 }
214 if(hitCount != theHighEnergyGenerator->GetWoundedNucleus()->GetMassNumber() )
215 {
216 theTransport->SetPrimaryProjectile(thePrimary);
217 theTransportResult =
218 theTransport->Propagate(theInitialResult,
219 theHighEnergyGenerator->GetWoundedNucleus());
220 if ( !theTransportResult ) {
221 G4cout << "G4TheoFSGenerator: null ptr from transport propagate " << G4endl;
222 throw G4HadronicException(__FILE__, __LINE__, "Null ptr from transport propagate");
223 }
224 }
225 else
226 {
227 theTransportResult = theDecay.Propagate(theInitialResult,
228 theHighEnergyGenerator->GetWoundedNucleus());
229 if ( theTransportResult == nullptr ) {
230 G4cout << "G4TheoFSGenerator: null ptr from decay propagate " << G4endl;
231 throw G4HadronicException(__FILE__, __LINE__, "Null ptr from decay propagate");
232 }
233 }
234 }
235 else
236 {
237 theTransport->SetPrimaryProjectile(thePrimary);
238 theTransportResult = theTransport->PropagateNuclNucl(theInitialResult,
239 theHighEnergyGenerator->GetWoundedNucleus(),
240 theProjectileNucleus);
241 if ( !theTransportResult ) {
242 G4cout << "G4TheoFSGenerator: null ptr from transport propagate " << G4endl;
243 throw G4HadronicException(__FILE__, __LINE__, "Null ptr from transport propagate");
244 }
245 }
246
247 // If enabled, apply the Cosmic Rays (CR) coalescence to the list of secondaries produced so far.
248 // This algorithm can form deuterons and antideuterons by coalescence of, respectively,
249 // proton-neutron and antiproton-antineutron pairs close in momentum space.
250 // This can be useful in particular for Cosmic Ray applications.
251 if ( G4HadronicParameters::Instance()->EnableCRCoalescence() ) {
252 if(nullptr == theCosmicCoalescence) {
253 theCosmicCoalescence = (G4CRCoalescence*)
254 G4HadronicInteractionRegistry::Instance()->FindModel("G4CRCoalescence");
255 if(nullptr == theCosmicCoalescence) {
256 theCosmicCoalescence = new G4CRCoalescence();
257 }
258 }
259 theCosmicCoalescence->SetP0Coalescence( thePrimary, theHighEnergyGenerator->GetModelName() );
260 theCosmicCoalescence->GenerateDeuterons( theTransportResult );
261 }
262
263 // Fill particle change
264 for(auto & ptr : *theTransportResult)
265 {
266 G4DynamicParticle * aNewDP =
267 new G4DynamicParticle(ptr->GetDefinition(),
268 ptr->GetTotalEnergy(),
269 ptr->GetMomentum());
270 G4HadSecondary aNew = G4HadSecondary(aNewDP);
271 G4double time = std::max(ptr->GetFormationTime(), 0.0);
272 aNew.SetTime(timePrimary + time);
273 aNew.SetCreatorModelID(ptr->GetCreatorModelID());
274 aNew.SetParentResonanceDef(ptr->GetParentResonanceDef());
275 aNew.SetParentResonanceID(ptr->GetParentResonanceID());
276 theParticleChange->AddSecondary(aNew);
277 delete ptr;
278 }
279
280 // some garbage collection
281 delete theTransportResult;
282
283 // Done
284 return theParticleChange;
285 }
286
287 std::pair<G4double, G4double> G4TheoFSGenerator::GetEnergyMomentumCheckLevels() const
288 {
289 if ( theHighEnergyGenerator ) {
290 return theHighEnergyGenerator->GetEnergyMomentumCheckLevels();
291 } else {
292 return std::pair<G4double, G4double>(DBL_MAX, DBL_MAX);
293 }
294 }
295