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