Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
G4AngularDistributionPP.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//
28#include "G4SystemOfUnits.hh"
29#include "Randomize.hh"
30#include "G4ios.hh"
31
32// Initialization of static data arrays:
34#include "Randomize.hh"
35
36
38{
39 G4int verboseLevel=1;
40
41 G4double ek= ((S - sqr(m_1) -sqr(m_2) )/(2*m_1) - m_1 )/GeV ; // kinetic energy in GeV
42
43 // Find energy bin
44
45 G4int je1 = 0;
46 G4int je2 = NENERGY - 1;
47 do {
48 G4int midBin = (je1 + je2)/2;
49 if (ek < elab[midBin])
50 je2 = midBin;
51 else
52 je1 = midBin;
53 } while (je2 - je1 > 1);
54 // G4int j;
55 //std::abs(ek-elab[je1]) < std::abs(ek-elab[je2]) ? j = je1 : j = je2;
56 G4double delab = elab[je2] - elab[je1];
57
58 // Sample the angle
59
60 G4float sample = G4UniformRand();
61 G4int ke1 = 0;
62 G4int ke2 = NANGLE - 1;
63 G4double dsig = sig[je2][0] - sig[je1][0];
64 G4double rc = dsig/delab;
65 G4double b = sig[je1][0] - rc*elab[je1];
66 G4double sigint1 = rc*ek + b;
67 G4double sigint2 = 0.;
68
69 if (verboseLevel > 1) G4cout << "sample=" << sample << G4endl
70 << ek << " " << ke1 << " " << ke2 << " "
71 << sigint1 << " " << sigint2 << G4endl;
72
73 do {
74 G4int midBin = (ke1 + ke2)/2;
75 dsig = sig[je2][midBin] - sig[je1][midBin];
76 rc = dsig/delab;
77 b = sig[je1][midBin] - rc*elab[je1];
78 G4double sigint = rc*ek + b;
79 if (sample < sigint) {
80 ke2 = midBin;
81 sigint2 = sigint;
82 }
83 else {
84 ke1 = midBin;
85 sigint1 = sigint;
86 }
87 if (verboseLevel > 1)G4cout << ke1 << " " << ke2 << " "
88 << sigint1 << " " << sigint2 << G4endl;
89 } while (ke2 - ke1 > 1);
90
91 // sigint1 and sigint2 should be recoverable from above loop
92
93 // G4double dsig = sig[je2][ke1] - sig[je1][ke1];
94 // G4double rc = dsig/delab;
95 // G4double b = sig[je1][ke1] - rc*elab[je1];
96 // G4double sigint1 = rc*ek + b;
97
98 // G4double dsig = sig[je2][ke2] - sig[je1][ke2];
99 // G4double rc = dsig/delab;
100 // G4double b = sig[je1][ke2] - rc*elab[je1];
101 // G4double sigint2 = rc*ek + b;
102
103 dsig = sigint2 - sigint1;
104 rc = 1./dsig;
105 b = ke1 - rc*sigint1;
106 G4double kint = rc*sample + b;
107 G4double theta = (0.5 + kint)*pi/180.;
108
109 // G4int k;
110 //std::abs(sample-sig[j][ke1]) < std::abs(sample-sig[j][ke2]) ? k = ke1 : k = ke2;
111 // G4double theta = (0.5 + k)*pi/180.;
112
113 if (verboseLevel > 1) {
114 G4cout << " energy bin " << je1 << " energy=" << elab[je1] << G4endl;
115 G4cout << " angle bin " << kint << " angle=" << theta/degree << G4endl;
116 }
117 G4double costh= std::cos(theta);
118// G4cout << "scattering angle hpw "<<costh<<" "<<s<<" "<<m1<<" "<<m2<<G4endl;
119 return costh;
120}
121
123{
124 return twopi * G4UniformRand();
125}
double G4double
Definition: G4Types.hh:64
float G4float
Definition: G4Types.hh:65
int G4int
Definition: G4Types.hh:66
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:53
virtual G4double Phi() const
virtual G4double CosTheta(G4double s, G4double m1, G4double m2) const
T sqr(const T &x)
Definition: templates.hh:145