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

Last change on this file since 1239 was 819, checked in by garnier, 16 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.