Geant4  10.00.p02
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 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
60 //
61 #include "G4EMDissociation.hh"
62 #include "G4PhysicalConstants.hh"
63 #include "G4SystemOfUnits.hh"
64 #include "G4Evaporation.hh"
65 #include "G4FermiBreakUp.hh"
66 #include "G4StatMF.hh"
67 #include "G4ParticleDefinition.hh"
68 #include "G4LorentzVector.hh"
69 #include "G4PhysicsFreeVector.hh"
71 #include "G4Proton.hh"
72 #include "G4Neutron.hh"
73 #include "G4ParticleTable.hh"
74 #include "G4IonTable.hh"
76 #include "G4DecayProducts.hh"
77 #include "G4DynamicParticle.hh"
78 #include "G4Fragment.hh"
80 #include "Randomize.hh"
81 #include "globals.hh"
82 
84 
85  // Send message to stdout to advise that the G4EMDissociation model is being
86  // used.
88 
89  // No de-excitation handler has been supplied - define the default handler.
91  G4Evaporation* theEvaporation = new G4Evaporation;
92  G4FermiBreakUp* theFermiBreakUp = new G4FermiBreakUp;
93  G4StatMF* theMF = new G4StatMF;
94  theExcitationHandler->SetEvaporation(theEvaporation);
95  theExcitationHandler->SetFermiModel(theFermiBreakUp);
100 
101  // This EM dissociation model needs access to the cross-sections held in
102  // G4EMDissociationCrossSection.
105 
106  // Set the minimum and maximum range for the model (despite nomanclature, this
107  // is in energy per nucleon number).
108  SetMinEnergy(100.0*MeV);
109  SetMaxEnergy(500.0*GeV);
110 
111  // Set the default verbose level to 0 - no output.
112  verboseLevel = 0;
113 }
114 
115 /*
116 G4EMDissociation::G4EMDissociation(const G4EMDissociation& emd)
117  : G4HadronicInteraction(emd)
118 {
119  if (emd.theExcitationHandler != 0) {
120  theExcitationHandler = new G4ExcitationHandler;
121  *theExcitationHandler = *emd.theExcitationHandler;
122  }
123 
124  handlerDefinedInternally = emd.handlerDefinedInternally;
125 
126  if (emd.dissociationCrossSection != 0) {
127  dissociationCrossSection = new G4EMDissociationCrossSection;
128  *dissociationCrossSection = *emd.dissociationCrossSection;
129  }
130 
131  if (emd.thePhotonSpectrum !- 0) {
132  thePhotonSpectrum = new G4EMDissociationSpectrum;
133  *thePhotonSpectrum = *emd.thePhotonSpectrum;
134 }
135 */
136 
138 {
139  // Send message to stdout to advise that the G4EMDissociation model is being
140  // used.
142 
143  theExcitationHandler = aExcitationHandler;
144  handlerDefinedInternally = false;
145 
146  // This EM dissociation model needs access to the cross-sections held in
147  // G4EMDissociationCrossSection.
150 
151  // Set the minimum and maximum range for the model (despite nomanclature, this
152  // is in energy per nucleon number)
153  SetMinEnergy(100.0*MeV);
154  SetMaxEnergy(500.0*GeV);
155  verboseLevel = 0;
156 }
157 
158 
161  // delete dissociationCrossSection;
162  // Cross section deleted by G4CrossSectionRegistry; don't do it here
163  // Bug reported by Gong Ding in Bug Report #1339
164  delete thePhotonSpectrum;
165 }
166 
167 
169  (const G4HadProjectile &theTrack, G4Nucleus &theTarget)
170 {
171  // The secondaries will be returned in G4HadFinalState &theParticleChange -
172  // initialise this.
173 
174  theParticleChange.Clear();
175  theParticleChange.SetStatusChange(stopAndKill);
176 
177  // Get relevant information about the projectile and target (A, Z) and
178  // energy/nuc, momentum, velocity, Lorentz factor and rest-mass of the
179  // projectile.
180 
181  const G4ParticleDefinition *definitionP = theTrack.GetDefinition();
182  const G4double AP = definitionP->GetBaryonNumber();
183  const G4double ZP = definitionP->GetPDGCharge();
184  G4LorentzVector pP = theTrack.Get4Momentum();
185  G4double E = theTrack.GetKineticEnergy()/AP;
186  G4double MP = theTrack.GetTotalEnergy() - E*AP;
187  G4double b = pP.beta();
188  G4double AT = theTarget.GetA_asInt();
189  G4double ZT = theTarget.GetZ_asInt();
191 
192  // Depending upon the verbosity level, output the initial information on the
193  // projectile and target
194  if (verboseLevel >= 2) {
195  G4cout.precision(6);
196  G4cout <<"########################################"
197  <<"########################################"
198  <<G4endl;
199  G4cout <<"IN G4EMDissociation" <<G4endl;
200  G4cout <<"Initial projectile A=" <<AP
201  <<", Z=" <<ZP
202  <<G4endl;
203  G4cout <<"Initial target A=" <<AT
204  <<", Z=" <<ZT
205  <<G4endl;
206  G4cout <<"Projectile momentum and Energy/nuc = " <<pP <<" ," <<E <<G4endl;
207  }
208 
209  // Initialise the variables which will be used with the phase-space decay and
210  // to boost the secondaries from the interaction.
211 
212  G4ParticleDefinition *typeNucleon = NULL;
213  G4ParticleDefinition *typeDaughter = NULL;
214  G4double Eg = 0.0;
215  G4double mass = 0.0;
216  G4ThreeVector boost = G4ThreeVector(0.0, 0.0, 0.0);
217 
218  // Determine the cross-sections at the giant dipole and giant quadrupole
219  // resonance energies for the projectile and then target. The information is
220  // initially provided in the G4PhysicsFreeVector individually for the E1
221  // and E2 fields. These are then summed.
222 
223  G4double bmin = thePhotonSpectrum->GetClosestApproach(AP, ZP, AT, ZT, b);
224  G4PhysicsFreeVector *crossSectionP = dissociationCrossSection->
225  GetCrossSectionForProjectile(AP, ZP, AT, ZT, b, bmin);
226  G4PhysicsFreeVector *crossSectionT = dissociationCrossSection->
227  GetCrossSectionForTarget(AP, ZP, AT, ZT, b, bmin);
228 
229  G4double totCrossSectionP = (*crossSectionP)[0]+(*crossSectionP)[1];
230  G4double totCrossSectionT = (*crossSectionT)[0]+(*crossSectionT)[1];
231 
232  // Now sample whether the interaction involved EM dissociation of the projectile
233  // or the target.
234 
235  if (G4UniformRand() <
236  totCrossSectionP / (totCrossSectionP + totCrossSectionT)) {
237 
238  // It was the projectile which underwent EM dissociation. Define the Lorentz
239  // boost to be applied to the secondaries, and sample whether a proton or a
240  // neutron was ejected. Then determine the energy of the virtual gamma ray
241  // which passed from the target nucleus ... this will be used to define the
242  // excitation of the projectile.
243 
244  mass = MP;
245  if (G4UniformRand() < dissociationCrossSection->
246  GetWilsonProbabilityForProtonDissociation (AP, ZP))
247  {
248  if (verboseLevel >= 2)
249  G4cout <<"Projectile underwent EM dissociation producing a proton"
250  <<G4endl;
251  typeNucleon = G4Proton::ProtonDefinition();
252  typeDaughter = G4ParticleTable::GetParticleTable()->
253  GetIon((G4int) ZP-1, (G4int) AP-1, 0.0);
254  }
255  else
256  {
257  if (verboseLevel >= 2)
258  G4cout <<"Projectile underwent EM dissociation producing a neutron"
259  <<G4endl;
260  typeNucleon = G4Neutron::NeutronDefinition();
261  typeDaughter = G4ParticleTable::GetParticleTable()->
262  GetIon((G4int) ZP, (G4int) AP-1, 0.0);
263  }
264  if (G4UniformRand() < (*crossSectionP)[0]/totCrossSectionP)
265  {
266  Eg = crossSectionP->GetLowEdgeEnergy(0);
267  if (verboseLevel >= 2)
268  G4cout <<"Transition type was E1" <<G4endl;
269  }
270  else
271  {
272  Eg = crossSectionP->GetLowEdgeEnergy(1);
273  if (verboseLevel >= 2)
274  G4cout <<"Transition type was E2" <<G4endl;
275  }
276 
277  // We need to define a Lorentz vector with the original momentum, but total
278  // energy includes the projectile and virtual gamma. This is then used
279  // to calculate the boost required for the secondaries.
280 
281  pP.setE(pP.e()+Eg);
282  boost = pP.findBoostToCM();
283  }
284  else
285  {
286  // It was the target which underwent EM dissociation. Sample whether a
287  // proton or a neutron was ejected. Then determine the energy of the virtual
288  // gamma ray which passed from the projectile nucleus ... this will be used to
289  // define the excitation of the target.
290 
291  mass = MT;
292  if (G4UniformRand() < dissociationCrossSection->
293  GetWilsonProbabilityForProtonDissociation (AT, ZT))
294  {
295  if (verboseLevel >= 2)
296  G4cout <<"Target underwent EM dissociation producing a proton"
297  <<G4endl;
298  typeNucleon = G4Proton::ProtonDefinition();
299  typeDaughter = G4ParticleTable::GetParticleTable()->
300  GetIon((G4int) ZT-1, (G4int) AT-1, 0.0);
301  }
302  else
303  {
304  if (verboseLevel >= 2)
305  G4cout <<"Target underwent EM dissociation producing a neutron"
306  <<G4endl;
307  typeNucleon = G4Neutron::NeutronDefinition();
308  typeDaughter = G4ParticleTable::GetParticleTable()->
309  GetIon((G4int) ZT, (G4int) AT-1, 0.0);
310  }
311  if (G4UniformRand() < (*crossSectionT)[0]/totCrossSectionT)
312  {
313  Eg = crossSectionT->GetLowEdgeEnergy(0);
314  if (verboseLevel >= 2)
315  G4cout <<"Transition type was E1" <<G4endl;
316  }
317  else
318  {
319  Eg = crossSectionT->GetLowEdgeEnergy(1);
320  if (verboseLevel >= 2)
321  G4cout <<"Transition type was E2" <<G4endl;
322  }
323 
324  // Add the projectile to theParticleChange, less the energy of the
325  // not-so-virtual gamma-ray. Not that at the moment, no lateral momentum
326  // is transferred between the projectile and target nuclei.
327 
328  G4ThreeVector v = pP.vect();
329  v.setMag(1.0);
330  G4DynamicParticle *changedP = new G4DynamicParticle
331  (const_cast<G4ParticleDefinition*>(definitionP), v, E*AP-Eg);
332  theParticleChange.AddSecondary (changedP);
333  if (verboseLevel >= 2)
334  {
335  G4cout <<"Projectile change:" <<G4endl;
336  changedP->DumpInfo();
337  }
338  }
339 
340  // Perform a two-body decay based on the restmass energy of the parent and
341  // gamma-ray, and the masses of the daughters. In the frame of reference of
342  // the nucles, the angular distribution is sampled isotropically, but the
343  // the nucleon and secondary nucleus are boosted if they've come from the
344  // projectile.
345 
346  G4double e = mass + Eg;
347  G4double mass1 = typeNucleon->GetPDGMass();
348  G4double mass2 = typeDaughter->GetPDGMass();
349  G4double pp = (e+mass1+mass2)*(e+mass1-mass2)*
350  (e-mass1+mass2)*(e-mass1-mass2)/(4.0*e*e);
351  if (pp < 0.0) {
352  pp = 1.0*eV;
353 // if (verboseLevel >`= 1)
354 // {
355 // G4cout <<"IN G4EMDissociation::ApplyYoursef" <<G4endl;
356 // G4cout <<"Error in mass of secondaries compared with primary:" <<G4endl;
357 // G4cout <<"Rest mass of primary = " <<mass <<" MeV" <<G4endl;
358 // G4cout <<"Virtual gamma energy = " <<Eg <<" MeV" <<G4endl;
359 // G4cout <<"Rest mass of secondary #1 = " <<mass1 <<" MeV" <<G4endl;
360 // G4cout <<"Rest mass of secondary #2 = " <<mass2 <<" MeV" <<G4endl;
361 // }
362  }
363  else
364  pp = std::sqrt(pp);
365  G4double costheta = 2.*G4UniformRand()-1.0;
366  G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
367  G4double phi = 2.0*pi*G4UniformRand()*rad;
368  G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
369  G4DynamicParticle *dynamicNucleon =
370  new G4DynamicParticle(typeNucleon, direction*pp);
371  dynamicNucleon->Set4Momentum(dynamicNucleon->Get4Momentum().boost(-boost));
372  G4DynamicParticle *dynamicDaughter =
373  new G4DynamicParticle(typeDaughter, -direction*pp);
374  dynamicDaughter->Set4Momentum(dynamicDaughter->Get4Momentum().boost(-boost));
375 
376  // The "decay" products have to be transferred to the G4HadFinalState object.
377  // Furthermore, the residual nucleus should be de-excited.
378 
379  theParticleChange.AddSecondary (dynamicNucleon);
380  if (verboseLevel >= 2) {
381  G4cout <<"Nucleon from the EMD process:" <<G4endl;
382  dynamicNucleon->DumpInfo();
383  }
384 
385  G4Fragment* theFragment = new
386  G4Fragment((G4int) typeDaughter->GetBaryonNumber(),
387  (G4int) typeDaughter->GetPDGCharge(), dynamicDaughter->Get4Momentum());
388 
389  if (verboseLevel >= 2) {
390  G4cout <<"Dynamic properties of the prefragment:" <<G4endl;
391  G4cout.precision(6);
392  dynamicDaughter->DumpInfo();
393  G4cout <<"Nuclear properties of the prefragment:" <<G4endl;
394  G4cout <<theFragment <<G4endl;
395  }
396 
397  G4ReactionProductVector* products =
398  theExcitationHandler->BreakItUp(*theFragment);
399  delete theFragment;
400  theFragment = NULL;
401 
402  G4DynamicParticle* secondary = 0;
403  G4ReactionProductVector::iterator iter;
404  for (iter = products->begin(); iter != products->end(); ++iter) {
405  secondary = new G4DynamicParticle((*iter)->GetDefinition(),
406  (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
407  theParticleChange.AddSecondary (secondary);
408  }
409  delete products;
410 
411  if (verboseLevel >= 2)
412  G4cout <<"########################################"
413  <<"########################################"
414  <<G4endl;
415 
416  return &theParticleChange;
417 }
418 
419 
421 {
422  G4cout <<G4endl;
423  G4cout <<" ****************************************************************"
424  <<G4endl;
425  G4cout <<" EM dissociation model for nuclear-nuclear interactions activated"
426  <<G4endl;
427  G4cout <<" (Written by QinetiQ Ltd for the European Space Agency)"
428  <<G4endl;
429  G4cout <<" ****************************************************************"
430  <<G4endl;
431  G4cout << G4endl;
432 
433  return;
434 }
435 
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4bool handlerDefinedInternally
static const double MeV
Definition: G4SIunits.hh:193
static G4double GetNuclearMass(const G4double A, const G4double Z)
CLHEP::Hep3Vector G4ThreeVector
const G4double pi
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
void DumpInfo(G4int mode=0) const
void SetMinEForMultiFrag(G4double anE)
int G4int
Definition: G4Types.hh:78
G4ExcitationHandler * theExcitationHandler
std::vector< G4ReactionProduct * > G4ReactionProductVector
void SetMinEnergy(G4double anEnergy)
G4EMDissociationCrossSection * dissociationCrossSection
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
void SetFermiModel(G4VFermiBreakUp *ptr)
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &, G4Nucleus &)
G4ErrorTarget * theTarget
Definition: errprop.cc:59
void SetMultiFragmentation(G4VMultiFragmentation *ptr)
static const double GeV
Definition: G4SIunits.hh:196
const G4LorentzVector & Get4Momentum() const
G4LorentzVector Get4Momentum() const
static const double rad
Definition: G4SIunits.hh:130
void Set4Momentum(const G4LorentzVector &momentum)
static const double eV
Definition: G4SIunits.hh:194
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
void SetMaxAandZForFermiBreakUp(G4int anA, G4int aZ)
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4EMDissociationSpectrum * thePhotonSpectrum
void SetEvaporation(G4VEvaporation *ptr)
void SetMaxEnergy(const G4double anEnergy)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
static G4Neutron * NeutronDefinition()
Definition: G4Neutron.cc:99
G4double GetTotalEnergy() const
CLHEP::HepLorentzVector G4LorentzVector