ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GammaTransition.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4GammaTransition.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 // * 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 //
27 // -------------------------------------------------------------------
28 // GEANT 4 class file
29 //
30 // CERN, Geneva, Switzerland
31 //
32 // File name: G4GammaTransition
33 //
34 // Author V.Ivanchenko 27 February 2015
35 //
36 // -------------------------------------------------------------------
37 
38 #include "G4GammaTransition.hh"
39 #include "G4AtomicShells.hh"
40 #include "Randomize.hh"
41 #include "G4RandomDirection.hh"
42 #include "G4Gamma.hh"
43 #include "G4Electron.hh"
44 #include "G4LorentzVector.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4PhysicalConstants.hh"
47 
49  : polarFlag(false), fDirection(0.,0.,0.), fTwoJMAX(10), fVerbose(0)
50 {}
51 
53 {}
54 
55 G4Fragment*
57  G4double newExcEnergy,
58  G4double mpRatio,
59  G4int JP1,
60  G4int JP2,
61  G4int MP,
62  G4int shell,
63  G4bool isDiscrete,
64  G4bool isGamma)
65 {
66  G4Fragment* result = nullptr;
67  G4double bond_energy = 0.0;
68 
69  if (!isGamma) {
70  if(0 <= shell) {
71  G4int Z = nucleus->GetZ_asInt();
72  if(Z <= 100) {
73  G4int idx = (G4int)shell;
75  bond_energy = G4AtomicShells::GetBindingEnergy(Z, idx);
76  }
77  }
78  }
79  G4double etrans = nucleus->GetExcitationEnergy() - newExcEnergy
80  - bond_energy;
81  if(fVerbose > 2) {
82  G4cout << "G4GammaTransition::GenerateGamma - Etrans(MeV)= "
83  << etrans << " Eexnew= " << newExcEnergy
84  << " Ebond= " << bond_energy << G4endl;
85  }
86  if(etrans <= 0.0) {
87  etrans += bond_energy;
88  bond_energy = 0.0;
89  }
90 
91  // Do complete Lorentz computation
92  G4LorentzVector lv = nucleus->GetMomentum();
93  G4double mass = nucleus->GetGroundStateMass() + newExcEnergy;
94 
95  // select secondary
97 
98  if(isGamma) { part = G4Gamma::Gamma(); }
99  else {
100  part = G4Electron::Electron();
101  G4int ne = std::max(nucleus->GetNumberOfElectrons() - 1, 0);
102  nucleus->SetNumberOfElectrons(ne);
103  }
104 
105  if(polarFlag && isDiscrete && JP1 <= fTwoJMAX) {
106  SampleDirection(nucleus, mpRatio, JP1, JP2, MP);
107  } else {
109  }
110 
111  G4double emass = part->GetPDGMass();
112 
113  // 2-body decay in rest frame
114  G4double ecm = lv.mag();
115  G4ThreeVector bst = lv.boostVector();
116  if(!isGamma) { ecm += (CLHEP::electron_mass_c2 - bond_energy); }
117 
118  //G4cout << "Ecm= " << ecm << " mass= " << mass << " emass= " << emass << G4endl;
119 
120  ecm = std::max(ecm, mass + emass);
121  G4double energy = 0.5*((ecm - mass)*(ecm + mass) + emass*emass)/ecm;
122  G4double mom = (emass > 0.0) ? std::sqrt((energy - emass)*(energy + emass))
123  : energy;
124 
125  // emitted gamma or e-
126  G4LorentzVector res4mom(mom * fDirection.x(),
127  mom * fDirection.y(),
128  mom * fDirection.z(), energy);
129  // residual
130  energy = std::max(ecm - energy, mass);
131  lv.set(-mom*fDirection.x(), -mom*fDirection.y(), -mom*fDirection.z(), energy);
132 
133  // Lab system transform for short lived level
134  lv.boost(bst);
135 
136  // modified primary fragment
137  nucleus->SetExcEnergyAndMomentum(newExcEnergy, lv);
138 
139  // gamma or e- are produced
140  res4mom.boost(bst);
141  result = new G4Fragment(res4mom, part);
142 
143  //G4cout << " DeltaE= " << e0 - lv.e() - res4mom.e() + emass
144  // << " Emass= " << emass << G4endl;
145  if(fVerbose > 2) {
146  G4cout << "G4GammaTransition::SampleTransition : " << *result << G4endl;
147  G4cout << " Left nucleus: " << *nucleus << G4endl;
148  }
149  return result;
150 }
151 
153  G4int twoJ1, G4int twoJ2, G4int mp)
154 {
155  G4double cosTheta, phi;
157  if(fVerbose > 2) {
158  G4cout << "G4GammaTransition::SampleDirection : 2J1= " << twoJ1
159  << " 2J2= " << twoJ2 << " ratio= " << ratio
160  << " mp= " << mp << G4endl;
161  G4cout << " Nucleus: " << *nuc << G4endl;
162  }
163  if(nullptr == np) {
164  cosTheta = 2*G4UniformRand() - 1.0;
165  phi = CLHEP::twopi*G4UniformRand();
166 
167  } else {
168  // PhotonEvaporation dataset:
169  // The multipolarity number with 1,2,3,4,5,6,7 representing E0,E1,M1,E2,M2,E3,M3
170  // monopole transition and 100*Nx+Ny representing multipolarity transition with
171  // Ny and Ny taking the value 1,2,3,4,5,6,7 referring to E0,E1,M1,E2,M2,E3,M3,..
172  // For example a M1+E2 transition would be written 304.
173  // M1 is the primary transition (L) and E2 is the secondary (L')
174 
175  G4double mpRatio = ratio;
176 
177  G4int L0 = 0, Lp = 0;
178  if (mp > 99) {
179  L0 = mp/200;
180  Lp = (mp%100)/2;
181  } else {
182  L0 = mp/2;
183  Lp = 0;
184  mpRatio = 0.;
185  }
186  fPolTrans.SampleGammaTransition(np, twoJ1, twoJ2, L0, Lp, mpRatio, cosTheta, phi);
187  }
188 
189  G4double sinTheta = std::sqrt((1.-cosTheta)*(1.+cosTheta));
190  fDirection.set(sinTheta*std::cos(phi),sinTheta*std::sin(phi),cosTheta);
191  if(fVerbose > 3) {
192  G4cout << "G4GammaTransition::SampleDirection done: " << fDirection << G4endl;
193  if(np) { G4cout << *np << G4endl; }
194  }
195 }