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