Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4EMDissociation.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// * *
21// * Parts of this code which have been developed by QinetiQ Ltd *
22// * under contract to the European Space Agency (ESA) are the *
23// * intellectual property of ESA. Rights to use, copy, modify and *
24// * redistribute this software for general public use are granted *
25// * in compliance with any licensing, distribution and development *
26// * policy adopted by the Geant4 Collaboration. This code has been *
27// * written by QinetiQ Ltd for the European Space Agency, under ESA *
28// * contract 17191/03/NL/LvH (Aurora Programme). *
29// * *
30// * By using, copying, modifying or distributing the software (or *
31// * any work based on the software) you agree to acknowledge its *
32// * use in resulting scientific publications, and indicate your *
33// * acceptance of all terms of the Geant4 Software license. *
34// ********************************************************************
35//
36// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37//
38// MODULE: G4EMDissociation.cc
39//
40// Version: B.1
41// Date: 15/04/04
42// Author: P R Truscott
43// Organisation: QinetiQ Ltd, UK
44// Customer: ESA/ESTEC, NOORDWIJK
45// Contract: 17191/03/NL/LvH
46//
47// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48//
49// CHANGE HISTORY
50// --------------
51//
52// 17 October 2003, P R Truscott, QinetiQ Ltd, UK
53// Created.
54//
55// 15 March 2004, P R Truscott, QinetiQ Ltd, UK
56// Beta release
57//
58// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
59////////////////////////////////////////////////////////////////////////////////
60//
61#include "G4EMDissociation.hh"
63#include "G4SystemOfUnits.hh"
65#include "G4LorentzVector.hh"
68#include "G4Proton.hh"
69#include "G4Neutron.hh"
70#include "G4IonTable.hh"
71#include "G4DecayProducts.hh"
72#include "G4DynamicParticle.hh"
73#include "G4Fragment.hh"
75#include "Randomize.hh"
76#include "globals.hh"
78
80 G4HadronicInteraction("EMDissociation"),
81 secID_projectileDissociation(-1), secID_targetDissociation(-1)
82{
83 // Send message to stdout to advise that the G4EMDissociation model is being
84 // used.
85 PrintWelcomeMessage();
86
87 // No de-excitation handler has been supplied - define the default handler.
88 theExcitationHandler = new G4ExcitationHandler;
89 theExcitationHandler->SetMinEForMultiFrag(5.0*MeV);
90 handlerDefinedInternally = true;
91
92 // This EM dissociation model needs access to the cross-sections held in
93 // G4EMDissociationCrossSection.
94 dissociationCrossSection = new G4EMDissociationCrossSection;
95 thePhotonSpectrum = new G4EMDissociationSpectrum;
96
97 // Set the minimum and maximum range for the model (despite nomanclature, this
98 // is in energy per nucleon number).
99 SetMinEnergy(100.0*MeV);
100 SetMaxEnergy(500.0*GeV);
101
102 // Set the default verbose level to 0 - no output.
103 verboseLevel = 0;
104
105 // Creator model ID for the secondaries created by this model
106 secID_projectileDissociation = G4PhysicsModelCatalog::GetModelID( "model_projectile" + GetModelName() );
107 secID_targetDissociation = G4PhysicsModelCatalog::GetModelID( "model_target" + GetModelName() );
108}
109
111 G4HadronicInteraction("EMDissociation"),
112 secID_projectileDissociation(-1), secID_targetDissociation(-1)
113{
114 // Send message to stdout to advise that the G4EMDissociation model is being
115 // used.
116 PrintWelcomeMessage();
117
118 theExcitationHandler = aExcitationHandler;
119 handlerDefinedInternally = false;
120
121 // This EM dissociation model needs access to the cross-sections held in
122 // G4EMDissociationCrossSection.
123 dissociationCrossSection = new G4EMDissociationCrossSection;
124 thePhotonSpectrum = new G4EMDissociationSpectrum;
125
126 // Set the minimum and maximum range for the model (despite nomanclature, this
127 // is in energy per nucleon number)
128 SetMinEnergy(100.0*MeV);
129 SetMaxEnergy(500.0*GeV);
130 verboseLevel = 0;
131
132 // Creator model ID for the secondaries created by this model
133 secID_projectileDissociation = G4PhysicsModelCatalog::GetModelID( "model_projectile" + GetModelName() );
134 secID_targetDissociation = G4PhysicsModelCatalog::GetModelID( "model_target" + GetModelName() );
135}
136
137
139 if (handlerDefinedInternally) delete theExcitationHandler;
140 // delete dissociationCrossSection;
141 // Cross section deleted by G4CrossSectionRegistry; don't do it here
142 // Bug reported by Gong Ding in Bug Report #1339
143 delete thePhotonSpectrum;
144}
145
146
148 (const G4HadProjectile &theTrack, G4Nucleus &theTarget)
149{
150 // The secondaries will be returned in G4HadFinalState &theParticleChange -
151 // initialise this.
152
155
156 // Get relevant information about the projectile and target (A, Z) and
157 // energy/nuc, momentum, velocity, Lorentz factor and rest-mass of the
158 // projectile.
159
160 const G4ParticleDefinition *definitionP = theTrack.GetDefinition();
161 const G4double AP = definitionP->GetBaryonNumber();
162 const G4double ZP = definitionP->GetPDGCharge();
163 G4LorentzVector pP = theTrack.Get4Momentum();
164 G4double E = theTrack.GetKineticEnergy()/AP;
165 G4double MP = theTrack.GetTotalEnergy() - E*AP;
166 G4double b = pP.beta();
167 G4double AT = theTarget.GetA_asInt();
168 G4double ZT = theTarget.GetZ_asInt();
170
171 // Depending upon the verbosity level, output the initial information on the
172 // projectile and target
173 if (verboseLevel >= 2) {
174 G4cout.precision(6);
175 G4cout <<"########################################"
176 <<"########################################"
177 <<G4endl;
178 G4cout <<"IN G4EMDissociation" <<G4endl;
179 G4cout <<"Initial projectile A=" <<AP
180 <<", Z=" <<ZP
181 <<G4endl;
182 G4cout <<"Initial target A=" <<AT
183 <<", Z=" <<ZT
184 <<G4endl;
185 G4cout <<"Projectile momentum and Energy/nuc = " <<pP <<" ," <<E <<G4endl;
186 }
187
188 // Initialise the variables which will be used with the phase-space decay and
189 // to boost the secondaries from the interaction.
190
191 G4ParticleDefinition *typeNucleon = NULL;
192 G4ParticleDefinition *typeDaughter = NULL;
193 G4double Eg = 0.0;
194 G4double mass = 0.0;
195 G4ThreeVector boost = G4ThreeVector(0.0, 0.0, 0.0);
196
197 // Determine the cross-sections at the giant dipole and giant quadrupole
198 // resonance energies for the projectile and then target. The information is
199 // initially provided in the G4PhysicsFreeVector individually for the E1
200 // and E2 fields. These are then summed.
201
202 G4double bmin = thePhotonSpectrum->GetClosestApproach(AP, ZP, AT, ZT, b);
203 G4PhysicsFreeVector *crossSectionP = dissociationCrossSection->
204 GetCrossSectionForProjectile(AP, ZP, AT, ZT, b, bmin);
205 G4PhysicsFreeVector *crossSectionT = dissociationCrossSection->
206 GetCrossSectionForTarget(AP, ZP, AT, ZT, b, bmin);
207
208 G4double totCrossSectionP = (*crossSectionP)[0]+(*crossSectionP)[1];
209 G4double totCrossSectionT = (*crossSectionT)[0]+(*crossSectionT)[1];
210
211 // Now sample whether the interaction involved EM dissociation of the projectile
212 // or the target.
213
214 G4int secID = -1; // Creator model ID for the secondaries
215 if (G4UniformRand() <
216 totCrossSectionP / (totCrossSectionP + totCrossSectionT)) {
217
218 // It was the projectile which underwent EM dissociation. Define the Lorentz
219 // boost to be applied to the secondaries, and sample whether a proton or a
220 // neutron was ejected. Then determine the energy of the virtual gamma ray
221 // which passed from the target nucleus ... this will be used to define the
222 // excitation of the projectile.
223
224 secID = secID_projectileDissociation;
225 mass = MP;
226 if (G4UniformRand() < dissociationCrossSection->
227 GetWilsonProbabilityForProtonDissociation (AP, ZP))
228 {
229 if (verboseLevel >= 2)
230 G4cout <<"Projectile underwent EM dissociation producing a proton"
231 <<G4endl;
232 typeNucleon = G4Proton::ProtonDefinition();
233 typeDaughter = G4IonTable::GetIonTable()->
234 GetIon((G4int) ZP-1, (G4int) AP-1, 0.0);
235 }
236 else
237 {
238 if (verboseLevel >= 2)
239 G4cout <<"Projectile underwent EM dissociation producing a neutron"
240 <<G4endl;
241 typeNucleon = G4Neutron::NeutronDefinition();
242 typeDaughter = G4IonTable::GetIonTable()->
243 GetIon((G4int) ZP, (G4int) AP-1, 0.0);
244 }
245 if (G4UniformRand() < (*crossSectionP)[0]/totCrossSectionP)
246 {
247 Eg = crossSectionP->GetLowEdgeEnergy(0);
248 if (verboseLevel >= 2)
249 G4cout <<"Transition type was E1" <<G4endl;
250 }
251 else
252 {
253 Eg = crossSectionP->GetLowEdgeEnergy(1);
254 if (verboseLevel >= 2)
255 G4cout <<"Transition type was E2" <<G4endl;
256 }
257
258 // We need to define a Lorentz vector with the original momentum, but total
259 // energy includes the projectile and virtual gamma. This is then used
260 // to calculate the boost required for the secondaries.
261
262 pP.setE( std::sqrt( pP.vect().mag2() + (mass + Eg)*(mass + Eg) ) );
263 boost = pP.findBoostToCM();
264 }
265 else
266 {
267 // It was the target which underwent EM dissociation. Sample whether a
268 // proton or a neutron was ejected. Then determine the energy of the virtual
269 // gamma ray which passed from the projectile nucleus ... this will be used to
270 // define the excitation of the target.
271
272 secID = secID_targetDissociation;
273 mass = MT;
274 if (G4UniformRand() < dissociationCrossSection->
275 GetWilsonProbabilityForProtonDissociation (AT, ZT))
276 {
277 if (verboseLevel >= 2)
278 G4cout <<"Target underwent EM dissociation producing a proton"
279 <<G4endl;
280 typeNucleon = G4Proton::ProtonDefinition();
281 typeDaughter = G4IonTable::GetIonTable()->
282 GetIon((G4int) ZT-1, (G4int) AT-1, 0.0);
283 }
284 else
285 {
286 if (verboseLevel >= 2)
287 G4cout <<"Target underwent EM dissociation producing a neutron"
288 <<G4endl;
289 typeNucleon = G4Neutron::NeutronDefinition();
290 typeDaughter = G4IonTable::GetIonTable()->
291 GetIon((G4int) ZT, (G4int) AT-1, 0.0);
292 }
293 if (G4UniformRand() < (*crossSectionT)[0]/totCrossSectionT)
294 {
295 Eg = crossSectionT->GetLowEdgeEnergy(0);
296 if (verboseLevel >= 2)
297 G4cout <<"Transition type was E1" <<G4endl;
298 }
299 else
300 {
301 Eg = crossSectionT->GetLowEdgeEnergy(1);
302 if (verboseLevel >= 2)
303 G4cout <<"Transition type was E2" <<G4endl;
304 }
305
306 // Add the projectile to theParticleChange, less the energy of the
307 // not-so-virtual gamma-ray. Not that at the moment, no lateral momentum
308 // is transferred between the projectile and target nuclei.
309
310 G4ThreeVector v = pP.vect();
311 v.setMag(1.0);
312 G4DynamicParticle *changedP = new G4DynamicParticle (definitionP, v, E*AP-Eg);
313 theParticleChange.AddSecondary (changedP, secID);
314 if (verboseLevel >= 2)
315 {
316 G4cout <<"Projectile change:" <<G4endl;
317 changedP->DumpInfo();
318 }
319 }
320
321 // Perform a two-body decay based on the restmass energy of the parent and
322 // gamma-ray, and the masses of the daughters. In the frame of reference of
323 // the nucles, the angular distribution is sampled isotropically, but the
324 // the nucleon and secondary nucleus are boosted if they've come from the
325 // projectile.
326
327 G4double e = mass + Eg;
328 G4double mass1 = typeNucleon->GetPDGMass();
329 G4double mass2 = typeDaughter->GetPDGMass();
330 G4double pp = (e+mass1+mass2)*(e+mass1-mass2)*
331 (e-mass1+mass2)*(e-mass1-mass2)/(4.0*e*e);
332 if (pp < 0.0) {
333 pp = 1.0*eV;
334// if (verboseLevel >`= 1)
335// {
336// G4cout <<"IN G4EMDissociation::ApplyYoursef" <<G4endl;
337// G4cout <<"Error in mass of secondaries compared with primary:" <<G4endl;
338// G4cout <<"Rest mass of primary = " <<mass <<" MeV" <<G4endl;
339// G4cout <<"Virtual gamma energy = " <<Eg <<" MeV" <<G4endl;
340// G4cout <<"Rest mass of secondary #1 = " <<mass1 <<" MeV" <<G4endl;
341// G4cout <<"Rest mass of secondary #2 = " <<mass2 <<" MeV" <<G4endl;
342// }
343 }
344 else
345 pp = std::sqrt(pp);
346 G4double costheta = 2.*G4UniformRand()-1.0;
347 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
348 G4double phi = 2.0*pi*G4UniformRand()*rad;
349 G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
350 G4DynamicParticle *dynamicNucleon =
351 new G4DynamicParticle(typeNucleon, direction*pp);
352 dynamicNucleon->Set4Momentum(dynamicNucleon->Get4Momentum().boost(-boost));
353 G4DynamicParticle *dynamicDaughter =
354 new G4DynamicParticle(typeDaughter, -direction*pp);
355 dynamicDaughter->Set4Momentum(dynamicDaughter->Get4Momentum().boost(-boost));
356
357 // The "decay" products have to be transferred to the G4HadFinalState object.
358 // Furthermore, the residual nucleus should be de-excited.
359
360 theParticleChange.AddSecondary (dynamicNucleon, secID);
361 if (verboseLevel >= 2) {
362 G4cout <<"Nucleon from the EMD process:" <<G4endl;
363 dynamicNucleon->DumpInfo();
364 }
365
366 G4Fragment* theFragment = new
367 G4Fragment(typeDaughter->GetBaryonNumber(),
368 G4lrint(typeDaughter->GetPDGCharge()/CLHEP::eplus),
369 dynamicDaughter->Get4Momentum());
370
371 if (verboseLevel >= 2) {
372 G4cout <<"Dynamic properties of the prefragment:" <<G4endl;
373 G4cout.precision(6);
374 dynamicDaughter->DumpInfo();
375 G4cout <<"Nuclear properties of the prefragment:" <<G4endl;
376 G4cout <<theFragment <<G4endl;
377 }
378
379 G4ReactionProductVector* products =
380 theExcitationHandler->BreakItUp(*theFragment);
381 delete theFragment;
382 theFragment = NULL;
383
384 G4DynamicParticle* secondary = 0;
385 G4ReactionProductVector::iterator iter;
386 for (iter = products->begin(); iter != products->end(); ++iter) {
387 secondary = new G4DynamicParticle((*iter)->GetDefinition(),
388 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
389 theParticleChange.AddSecondary (secondary, secID);
390 }
391 delete products;
392
393 delete crossSectionP;
394 delete crossSectionT;
395
396 if (verboseLevel >= 2)
397 G4cout <<"########################################"
398 <<"########################################"
399 <<G4endl;
400
401 return &theParticleChange;
402}
403
404
405void G4EMDissociation::PrintWelcomeMessage ()
406{
407 G4cout <<G4endl;
408 G4cout <<" ****************************************************************"
409 <<G4endl;
410 G4cout <<" EM dissociation model for nuclear-nuclear interactions activated"
411 <<G4endl;
412 G4cout <<" (Written by QinetiQ Ltd for the European Space Agency)"
413 <<G4endl;
414 G4cout <<" ****************************************************************"
415 <<G4endl;
416 G4cout << G4endl;
417
418 return;
419}
420
@ stopAndKill
std::vector< G4ReactionProduct * > G4ReactionProductVector
CLHEP::Hep3Vector G4ThreeVector
double G4double
Definition: G4Types.hh:83
int G4int
Definition: G4Types.hh:85
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
double mag2() const
void setMag(double)
Definition: ThreeVector.cc:20
HepLorentzVector & boost(double, double, double)
Hep3Vector vect() const
Hep3Vector findBoostToCM() const
void DumpInfo(G4int mode=0) const
G4LorentzVector Get4Momentum() const
void Set4Momentum(const G4LorentzVector &momentum)
G4double GetClosestApproach(const G4double, const G4double, G4double, G4double, G4double)
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &, G4Nucleus &)
G4ReactionProductVector * BreakItUp(const G4Fragment &theInitialState)
void SetMinEForMultiFrag(G4double anE)
void SetStatusChange(G4HadFinalStateStatus aS)
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
const G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
const G4LorentzVector & Get4Momentum() const
G4double GetTotalEnergy() const
void SetMinEnergy(G4double anEnergy)
const G4String & GetModelName() const
void SetMaxEnergy(const G4double anEnergy)
static G4IonTable * GetIonTable()
Definition: G4IonTable.cc:170
static G4Neutron * NeutronDefinition()
Definition: G4Neutron.cc:98
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4int GetA_asInt() const
Definition: G4Nucleus.hh:99
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:105
G4double GetPDGCharge() const
static G4int GetModelID(const G4int modelIndex)
G4double GetLowEdgeEnergy(const std::size_t index) const
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:87
int G4lrint(double ad)
Definition: templates.hh:134