source: trunk/source/processes/hadronic/models/em_dissociation/src/G4EMDissociation.cc@ 1036

Last change on this file since 1036 was 819, checked in by garnier, 18 years ago

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