Geant4  10.02.p01
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 94676 2015-12-02 09:51:20Z gunter $
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  _tolerance = 0.1*CLHEP::keV;
81 }
82 
84 {
85  delete _transition;
86 }
87 
89  G4Fragment* nucleus)
90 {
91  if (_verbose > 1) { G4cout << "G4VGammaDeexcitation::DoChain" << G4endl; }
92 
93  if(CanDoTransition(nucleus)) {
94  for(size_t i=0; i<100; ++i) {
96  G4Fragment* gamma = GenerateGamma(nucleus);
97  if (gamma) { products->push_back(gamma); }
98  else { break; }
99  //G4cout << i << ". Egamma(MeV)= " << gamma->GetMomentum().e()
100  // << "; new Eex(MeV)= " << nucleus->GetExcitationEnergy()
101  // << G4endl;
102  if(nucleus->GetExcitationEnergy() <= _tolerance) { break; }
103  }
104  }
105  if (_verbose > 1) {
106  G4cout << "G4VGammaDeexcitation::DoChain - end" << G4endl;
107  }
108 }
109 
111 {
112  G4Fragment * thePhoton = 0;
113  _vSN = -1;
114 
115  _transition->SelectGamma(); // it can be conversion electron too
116 
117  G4double etrans = _transition->GetGammaEnergy();
118 
119  //L.Desorgher 05/01/2015 need to add the bond energy for correct
120  // computation of a transition in case of ICM
121  G4DiscreteGammaTransition* dtransition =
122  dynamic_cast <G4DiscreteGammaTransition*> (_transition);
123  G4double bond_energy=0.;
124 
125  if (dtransition && !dtransition->IsAGamma()) {
126  bond_energy = dtransition->GetBondEnergy();
127  }
128  etrans += bond_energy;
129  //G4cout << "G4VGammaDeexcitation::GenerateGamma - Etrans(MeV)= "
130  // << etrans << G4endl;
131  if(etrans <= 0.0) { return thePhoton; }
132 
133  // final excitation
134  G4double excitation = aNucleus->GetExcitationEnergy() - etrans;
135  if(excitation <= _tolerance) { excitation = 0.0; }
136 
138  if (_verbose > 1) {
139  G4cout << "G4VGammaDeexcitation::GenerateGamma - Edeexc(MeV)= "
140  << etrans << "; Time(ns)= " << gammaTime/CLHEP::ns
141  << "; left Eexc(MeV)= " << excitation << G4endl;
142  }
143 
144  // Do complete Lorentz computation
145  G4LorentzVector lv = aNucleus->GetMomentum();
146  G4double Mass = aNucleus->GetGroundStateMass() + excitation;
147 
148  // select secondary
150 
151  if (dtransition && !dtransition->IsAGamma() ) {
152  gamma = G4Electron::Electron();
153  _vSN = dtransition->GetOrbitNumber();
155  //L. Desorgher 05/01/2015 need to remove atomic bond energy
156  // of the IC electron
157  lv += G4LorentzVector(0.0,0.0,0.0,
158  CLHEP::electron_mass_c2 - bond_energy);
159  }
160 
161  G4double cosTheta = 1. - 2. * G4UniformRand();
162  G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
163  G4double phi = twopi * G4UniformRand();
164 
165  G4double eMass = gamma->GetPDGMass();
166  G4LorentzVector Gamma4P;
167  /*
168  G4cout << " Mass= " << eMass << " t= " << gammaTime
169  << " tlim= " << _timeLimit << G4endl;
170  */
171  // 2-body decay in rest frame
172  G4double Ecm = lv.mag();
173  G4ThreeVector bst = lv.boostVector();
174 
175  G4double GammaEnergy = 0.5*((Ecm - Mass)*(Ecm + Mass) + eMass*eMass)/Ecm;
176  if(GammaEnergy < eMass) { GammaEnergy = eMass; }
177 
178  G4double mom = std::sqrt((GammaEnergy - eMass)*(GammaEnergy + eMass));
179  Gamma4P.set(mom * sinTheta * std::cos(phi),
180  mom * sinTheta * std::sin(phi),
181  mom * cosTheta, GammaEnergy);
182 
183  // Lab system in normal case (_timeLimit = DBL_MAX)
184  if(gammaTime <= _timeLimit) {
185  Gamma4P.boost(bst);
186  lv -= Gamma4P;
187  } else {
188  // In exceptional case sample decay at rest at not correct position
189  // of stopping ion, 4-momentum balance is breaked but gamma energy
190  // is correct
191  lv -= Gamma4P;
192  G4double E = lv.e();
193  G4double P2= (E - Mass)*(E + Mass);
194  G4ThreeVector v = lv.vect().unit();
195  G4double p = 0.0;
196  if(P2 > 0.0) { p = std::sqrt(P2); }
197  else { E = Mass; }
198  lv.set(v.x()*p, v.y()*p, v.z()*p, E);
199  }
200 
201  // modified primary fragment
202  gammaTime += aNucleus->GetCreationTime();
203  aNucleus->SetMomentum(lv);
204  aNucleus->SetCreationTime(gammaTime);
205 
206  // gamma or e- are produced
207  thePhoton = new G4Fragment(Gamma4P,gamma);
208  thePhoton->SetCreationTime(gammaTime);
209 
210  //G4cout << "G4VGammaDeexcitation::GenerateGamma : " << thePhoton << G4endl;
211  //G4cout << " Left nucleus: " << aNucleus << G4endl;
212  return thePhoton;
213 }
214 
virtual G4double GetGammaEnergy()=0
CLHEP::Hep3Vector G4ThreeVector
G4ElectronOccupancy _electronO
virtual void SetEnergyFrom(G4double energy)=0
virtual G4bool CanDoTransition(G4Fragment *aNucleus)=0
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
virtual void SelectGamma()=0
G4double GetCreationTime() const
Definition: G4Fragment.hh:405
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:284
void SetMomentum(const G4LorentzVector &value)
Definition: G4Fragment.hh:289
static const double twopi
Definition: G4SIunits.hh:75
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:63
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:273
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4Fragment * GenerateGamma(G4Fragment *nucleus)
void SetCreationTime(G4double time)
Definition: G4Fragment.hh:410
G4double GetPDGMass() const
G4VGammaTransition * _transition
static const G4double * P2[nN]
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
static const double keV
Definition: G4SIunits.hh:213
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
#define ns
Definition: xmlparse.cc:614
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:268
G4int RemoveElectron(G4int orbit, G4int number=1)
void DoChain(G4FragmentVector *, G4Fragment *nucleus)
CLHEP::HepLorentzVector G4LorentzVector
virtual G4double GetGammaCreationTime()=0