ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DeltaAngle.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DeltaAngle.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 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4DeltaAngle
33 //
34 // Author: Vladimir Ivantcheko
35 //
36 // Creation date: 23 August 2013
37 //
38 // Modifications:
39 //
40 // Class Description:
41 //
42 // Delta-electron Angular Distribution Generation
43 //
44 // Class Description: End
45 //
46 // -------------------------------------------------------------------
47 //
48 
49 #include "G4DeltaAngle.hh"
50 #include "G4PhysicalConstants.hh"
51 #include "Randomize.hh"
52 #include "G4ParticleDefinition.hh"
53 #include "G4Electron.hh"
54 #include "G4AtomicShells.hh"
55 #include "G4SystemOfUnits.hh"
56 #include "G4Log.hh"
57 
58 using namespace std;
59 
61  : G4VEmAngularDistribution("deltaVI")
62 {
64  nprob = 26;
65  prob.resize(nprob,0.0);
66  fShellIdx = -1;
67 }
68 
70 {}
71 
74  G4double kinEnergyFinal, G4int Z, G4int idx,
75  const G4Material* mat)
76 {
77  fShellIdx = idx;
78  return SampleDirection(dp, kinEnergyFinal,Z, mat);
79 }
80 
83  G4double kinEnergyFinal, G4int Z,
84  const G4Material*)
85 {
88 
89  // if idx is not properly defined sample shell index
90  if(idx < 0 || idx >= nShells) {
91  if(nShells> nprob) {
92  nprob = nShells;
93  prob.resize(nprob,0.0);
94  }
95  G4double sum = 0.0;
96  for(idx=0; idx<nShells; ++idx) {
99  prob[idx] = sum;
100  }
101  sum *= G4UniformRand();
102  for(idx=0; idx<nShells; ++idx) {
103  if(sum <= prob[idx]) { break; }
104  }
105  }
107  G4double cost;
108  /*
109  G4cout << "E(keV)= " << kinEnergyFinal/keV
110  << " Ebind(keV)= " << bindingEnergy
111  << " idx= " << idx << " nShells= " << nShells << G4endl;
112  */
113  G4int n = 0;
114  G4bool isOK = false;
115  static const G4int nmax = 100;
116  do {
117  ++n;
118  // the atomic electron
120  G4double eKinEnergy = bindingEnergy*x;
121  G4double ePotEnergy = bindingEnergy*(1.0 + x);
122  G4double e = kinEnergyFinal + ePotEnergy + electron_mass_c2;
123  G4double p = sqrt((e + electron_mass_c2)*(e - electron_mass_c2));
124 
125  G4double totEnergy = dp->GetTotalEnergy();
126  G4double totMomentum = dp->GetTotalMomentum();
127  if(dp->GetParticleDefinition() == fElectron) {
128  totEnergy += ePotEnergy;
129  totMomentum = sqrt((totEnergy + electron_mass_c2)
130  *(totEnergy - electron_mass_c2));
131  }
132 
133  G4double eTotEnergy = eKinEnergy + electron_mass_c2;
134  G4double eTotMomentum = sqrt(eKinEnergy*(eTotEnergy + electron_mass_c2));
135  G4double costet = 2*G4UniformRand() - 1;
136  G4double sintet = sqrt((1 - costet)*(1 + costet));
137 
138  cost = 1.0;
139  if(n >= nmax) {
140  /*
141  G4ExceptionDescription ed;
142  ed << "### G4DeltaAngle Warning: " << n
143  << " iterations - stop the loop with cost= 1.0 "
144  << " for " << dp->GetDefinition()->GetParticleName() << "\n"
145  << " Ekin(MeV)= " << dp->GetKineticEnergy()/MeV
146  << " Efinal(MeV)= " << kinEnergyFinal/MeV
147  << " Ebinding(MeV)= " << bindingEnergy/MeV;
148  G4Exception("G4DeltaAngle::SampleDirection","em0044",
149  JustWarning, ed,"");
150  */
151  if(0.0 == bindingEnergy) { isOK = true; }
152  bindingEnergy = 0.0;
153  }
154 
155  G4double x0 = p*(totMomentum + eTotMomentum*costet);
156  /*
157  G4cout << " x0= " << x0 << " p= " << p
158  << " ptot= " << totMomentum << " pe= " << eTotMomentum
159  << " e= " << e << " totMom= " << totMomentum
160  << G4endl;
161  */
162  if(x0 > 0.0) {
163  G4double x1 = p*eTotMomentum*sintet;
164  G4double x2 = totEnergy*(eTotEnergy - e) - e*eTotEnergy
165  - totMomentum*eTotMomentum*costet + electron_mass_c2*electron_mass_c2;
166  G4double y = -x2/x0;
167  if(std::abs(y) <= 1.0) {
168  cost = -(x2 + x1*sqrt(1. - y*y))/x0;
169  if(std::abs(cost) <= 1.0) { isOK = true; }
170  else { cost = 1.0; }
171  }
172 
173  /*
174  G4cout << " Ekin(MeV)= " << dp->GetKineticEnergy()
175  << " e1(keV)= " << eKinEnergy/keV
176  << " e2(keV)= " << (e - electron_mass_c2)/keV
177  << " 1-cost= " << 1-cost
178  << " x0= " << x0 << " x1= " << x1 << " x2= " << x2
179  << G4endl;
180  */
181  }
182 
183  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
184  } while(!isOK);
185 
186  G4double sint = sqrt((1 - cost)*(1 + cost));
188 
189  fLocalDirection.set(sint*cos(phi), sint*sin(phi), cost);
191 
192  return fLocalDirection;
193 }
194 
196 {}