Geant4_10
G4VGammaDeexcitation.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 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4VGammaDeexcitation.cc 67983 2013-03-13 10:42:03Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 // GEANT 4 class file
30 //
31 // CERN, Geneva, Switzerland
32 //
33 // File name: G4VGammaDeexcitation
34 //
35 // Author: Maria Grazia Pia (pia@genova.infn.it)
36 //
37 // Creation date: 23 October 1998
38 //
39 // Modifications:
40 //
41 // 21 Nov 2001, Fan Lei (flei@space.qinetiq.com)
42 // Modified GenerateGamma() and UpdateUncleus() for implementation
43 // of Internal Conversion processs
44 //
45 // 15 April 1999, Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
46 // Added creation time evaluation for products of evaporation
47 //
48 // 19 April 2010, J. M. Quesada calculations in CM system
49 // pending final boost to lab system (not critical)
50 //
51 // 23 April 2010, V.Ivanchenko rewite kinematic part using PDG formula
52 // for 2-body decay
53 //
54 // 07 May 2011, V.Ivanchenko implement check ICM flag - produce or not e-
55 //
56 // -------------------------------------------------------------------
57 
58 #include "G4VGammaDeexcitation.hh"
59 
60 #include "globals.hh"
61 #include "G4PhysicalConstants.hh"
62 #include "Randomize.hh"
63 #include "G4Gamma.hh"
64 #include "G4Electron.hh"
65 #include "G4LorentzVector.hh"
66 #include "G4VGammaTransition.hh"
67 #include "G4Fragment.hh"
68 #include "G4FragmentVector.hh"
69 
70 #include "G4ParticleTable.hh"
71 #include "G4IonTable.hh"
72 
74 
75 
76 G4VGammaDeexcitation::G4VGammaDeexcitation(): _transition(0), _verbose(0),
77  _electronO (0), _vSN(-1)
78 {
79  _nucleus = 0;
80  fTimeLimit = DBL_MAX;
81 }
82 
84 {
85  if (_transition != 0) { delete _transition; }
86 }
87 
89 {
90  Initialize();
91  G4FragmentVector* products = new G4FragmentVector();
92 
93  if (CanDoTransition())
94  {
95  G4Fragment* gamma = GenerateGamma();
96  if (gamma != 0) { products->push_back(gamma); }
97  }
98 
99  if (_verbose > 1) {
100  G4cout << "G4VGammaDeexcitation::DoTransition - Transition deleted " << G4endl;
101  }
102 
103  return products;
104 }
105 
107 {
108  if (_verbose > 1) { G4cout << "G4VGammaDeexcitation::DoChain" << G4endl; }
109  const G4double tolerance = CLHEP::keV;
110 
111  Initialize();
112  G4FragmentVector* products = new G4FragmentVector();
113 
114  while (CanDoTransition())
115  {
117  G4Fragment* gamma = GenerateGamma();
118  if (gamma != 0)
119  {
120  products->push_back(gamma);
121  //G4cout << "Eex(keV)= " << _nucleus->GetExcitationEnergy()/keV << G4endl;
122  if(_nucleus->GetExcitationEnergy() <= tolerance) { break; }
123  Update();
124  }
125  }
126 
127  if (_verbose > 1) {
128  G4cout << "G4VGammaDeexcitation::DoChain - Transition deleted, end of chain " << G4endl;
129  }
130 
131  return products;
132 }
133 
135 {
136  // 23/04/10 V.Ivanchenko rewrite complitely
137  G4double eGamma = 0.;
138 
139  if (_transition) {
140  _transition->SelectGamma(); // it can be conversion electron too
141  eGamma = _transition->GetGammaEnergy();
142  //G4cout << "G4VGammaDeexcitation::GenerateGamma - Egam(MeV)= "
143  // << eGamma << G4endl;
144  if(eGamma <= 0.0) { return 0; }
145  } else { return 0; }
146 
147  G4double excitation = _nucleus->GetExcitationEnergy() - eGamma;
148  if(excitation < 0.0) { excitation = 0.0; }
149  if (_verbose > 1)
150  {
151  G4cout << "G4VGammaDeexcitation::GenerateGamma - Edeexc(MeV)= " << eGamma
152  << " ** left Eexc(MeV)= " << excitation
153  << G4endl;
154  }
155 
157 
158  // Do complete Lorentz computation
159  G4LorentzVector lv = _nucleus->GetMomentum();
160  G4double Mass = _nucleus->GetGroundStateMass() + excitation;
161 
162  // select secondary
164 
165  G4DiscreteGammaTransition* dtransition =
166  dynamic_cast <G4DiscreteGammaTransition*> (_transition);
167 
168  G4bool eTransition = false;
169  if (dtransition && !( dtransition->IsAGamma()) ) {
170  eTransition = true;
171  gamma = G4Electron::Electron();
172  _vSN = dtransition->GetOrbitNumber();
173  _electronO.RemoveElectron(_vSN);
174  lv += G4LorentzVector(0.0,0.0,0.0,CLHEP::electron_mass_c2 - dtransition->GetBondEnergy());
175  }
176 
177  G4double cosTheta = 1. - 2. * G4UniformRand();
178  G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
179  G4double phi = twopi * G4UniformRand();
180 
181  G4double eMass = gamma->GetPDGMass();
182  G4LorentzVector Gamma4P;
183 
184  //G4cout << "Egamma= " << eGamma << " Mass= " << eMass << " t= " << gammaTime
185  // << " tlim= " << fTimeLimit << G4endl;
186 
187  if(gammaTime > fTimeLimit) {
188  // shortcut for long lived levels
189  // not correct position of stopping ion gamma emission
190  // 4-momentum balance is breaked
191  G4double e = eGamma + eMass;
192  G4double mom = std::sqrt(eGamma*(eGamma + 2*eMass));
193  Gamma4P.set(mom * sinTheta * std::cos(phi),
194  mom * sinTheta * std::sin(phi),
195  mom * cosTheta, e);
196  lv -= Gamma4P;
197  e = lv.e();
198  if(e < Mass) { e = Mass; }
199  mom = std::sqrt((e - Mass)*(e + Mass));
200  G4ThreeVector v = lv.vect().unit();
201  lv.set(mom*v.x(), mom*v.y(), mom*v.z(), e);
202 
203  } else {
204  // 2-body decay in rest frame
205  G4double Ecm = lv.mag();
206  G4ThreeVector bst = lv.boostVector();
207 
208  G4double GammaEnergy = 0.5*((Ecm - Mass)*(Ecm + Mass) + eMass*eMass)/Ecm;
209  if(GammaEnergy <= eMass) { return 0; }
210 
211  G4double mom = std::sqrt((GammaEnergy - eMass)*(GammaEnergy + eMass));
212  Gamma4P.set(mom * sinTheta * std::cos(phi),
213  mom * sinTheta * std::sin(phi),
214  mom * cosTheta,
215  GammaEnergy);
216 
217  Gamma4P.boost(bst);
218  lv -= Gamma4P;
219  }
220 
221  // modified primary fragment
222  gammaTime += _nucleus->GetCreationTime();
223 
224  _nucleus->SetMomentum(lv);
225  _nucleus->SetCreationTime(gammaTime);
226 
227  // e- is not produced
228  if(eTransition && !dtransition->GetICM()) { return 0; }
229 
230  // gamma or e- are produced
231  G4Fragment * thePhoton = new G4Fragment(Gamma4P,gamma);
232  thePhoton->SetCreationTime(gammaTime);
233 
234  //G4cout << "G4VGammaDeexcitation::GenerateGamma : " << thePhoton << G4endl;
235  //G4cout << " Left nucleus: " << _nucleus << G4endl;
236  return thePhoton;
237 }
238 
240 {
241  if (_transition != 0)
242  {
243  delete _transition;
244  _transition = 0;
245  if (_verbose > 1) {
246  G4cout << "G4VGammaDeexcitation::Update - Transition deleted " << G4endl;
247  }
248  }
249 
251  if (_transition != 0)
252  {
254  // if ( _vSN != -1) (dynamic_cast <G4DiscreteGammaTransition*> (_transition))->SetICM(false);
255  // the above line is commented out for bug fix #952. It was intruduced for reason that
256  // the k-shell electron is most likely one to be kicked out and there is no time for
257  // the atom to deexcite before the next IC. But this limitation is causing other problems as
258  // reported in #952
259  }
260 
261  return;
262 }
Hep3Vector boostVector() const
virtual G4bool CanDoTransition()=0
virtual G4double GetGammaEnergy()=0
double x() const
virtual void SetEnergyFrom(G4double energy)=0
double z() const
virtual G4VGammaTransition * CreateTransition()=0
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
virtual void SelectGamma()=0
G4double GetCreationTime() const
Definition: G4Fragment.hh:378
double mag() const
bool G4bool
Definition: G4Types.hh:79
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:251
HepLorentzVector & boost(double, double, double)
void SetMomentum(const G4LorentzVector &value)
Definition: G4Fragment.hh:256
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:64
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:240
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4FragmentVector * DoChain()
tuple v
Definition: test.py:18
void SetCreationTime(G4double time)
Definition: G4Fragment.hh:383
void set(double x, double y, double z, double t)
G4double GetPDGMass() const
G4FragmentVector * DoTransition()
Hep3Vector unit() const
double y() const
G4VGammaTransition * _transition
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:235
G4int RemoveElectron(G4int orbit, G4int number=1)
CLHEP::HepLorentzVector G4LorentzVector
virtual G4double GetGammaCreationTime()=0