Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4INCLEtaNElasticChannel.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// INCL++ intra-nuclear cascade model
27// Alain Boudard, CEA-Saclay, France
28// Joseph Cugnon, University of Liege, Belgium
29// Jean-Christophe David, CEA-Saclay, France
30// Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31// Sylvie Leray, CEA-Saclay, France
32// Davide Mancusi, CEA-Saclay, France
33//
34#define INCLXX_IN_GEANT4_MODE 1
35
36#include "globals.hh"
37
38#include "G4INCLEtaNElasticChannel.hh"
39#include "G4INCLKinematicsUtils.hh"
40#include "G4INCLBinaryCollisionAvatar.hh"
41#include "G4INCLRandom.hh"
42#include "G4INCLGlobals.hh"
43#include "G4INCLLogger.hh"
44
45namespace G4INCL {
46
47 EtaNElasticChannel::EtaNElasticChannel(Particle *p1, Particle *p2)
48 : particle1(p1), particle2(p2)
49 {
50
51 }
52
53 EtaNElasticChannel::~EtaNElasticChannel(){
54
55 }
56
57 void EtaNElasticChannel::fillFinalState(FinalState *fs) {
58 Particle * nucleon;
59 Particle * eta;
60 if(particle1->isNucleon()) {
61 nucleon = particle1;
62 eta = particle2;
63 } else {
64 nucleon = particle2;
65 eta = particle1;
66 }
67
68 G4double plab=KinematicsUtils::momentumInLab(particle1, particle2);
69
70 G4double sh=nucleon->getEnergy()+eta->getEnergy();
71 G4double mn=nucleon->getMass();
72 G4double me=eta->getMass();
73 G4double en=(sh*sh+mn*mn-me*me)/(2*sh);
74 nucleon->setEnergy(en);
75 G4double ee=std::sqrt(en*en-mn*mn+me*me);
76 eta->setEnergy(ee);
77 G4double pn=std::sqrt(en*en-mn*mn);
78
79 ThreeVector mom_nucleon;
80
81 if (plab < 250.) {
82// Isotropy
83 mom_nucleon = Random::normVector(pn);
84 }
85
86// From Kamano
87 else {
88
89 const G4double pi=std::acos(-1.0);
90 G4double x1;
91 G4double u1;
92 G4double fteta;
93 G4double teta;
94 G4double fi;
95
96 G4double a0;
97 G4double a1;
98 G4double a2;
99 G4double a3;
100 G4double a4;
101 G4double a5;
102 G4double a6;
103
104 if (plab > 1400.) plab=1400.; // no information on angular distributions above plab=1400 MeV
105 G4double p6=std::pow(plab, 6);
106 G4double p5=std::pow(plab, 5);
107 G4double p4=std::pow(plab, 4);
108 G4double p3=std::pow(plab, 3);
109 G4double p2=std::pow(plab, 2);
110 G4double p1=plab;
111
112// a6
113 if (plab < 300.) {
114 a6=-8.384000E-08*p1 - 1.15452E-04;
115 }
116 else if (plab < 500.){
117 a6=1.593966E-13*p4 - 2.619560E-10*p3 + 1.564701E-07*p2 - 3.986627E-05*p1 + 3.622575E-03;
118 }
119 else {
120 a6=6.143615E-20*p6 - 3.157181E-16*p5 + 6.348289E-13*p4 - 6.117961E-10*p3 + 2.764542E-07*p2 - 4.391048E-05*p1 - 1.443857E-03;
121 }
122// a5
123 if (plab < 650.) {
124 a5=-9.021076E-18*p6 + 2.176771E-14*p5 - 2.136095E-11*p4 + 1.100580E-08*p3 - 3.150857E-06*p2 + 4.761016E-04*p1 - 2.969608E-02;
125 }
126 else if (plab < 950.){
127 a5=4.424756E-18*p6 - 1.756295E-14*p5 + 2.625428E-11*p4 - 1.678272E-08*p3 + 2.227237E-06*p2 + 2.146666E-03*p1 - 7.065712E-01;
128 }
129 else {
130 a5=2.209585E-19*p6 - 1.546647E-15*p5 + 4.578142E-12*p4 - 7.303856E-09*p3 + 6.604074E-06*p2 - 3.205628E-03*p1 + 6.534893E-01;
131 }
132// a4
133 if (plab < 700.) {
134 a4=4.826684E-17*p6 - 1.534471E-13*p5 + 1.907868E-10*p4 - 1.192317E-07*p3 + 3.988902E-05*p2 - 6.822100E-03*p1 + 4.684685E-01;
135 }
136 else {
137 a4=-3.245143E-18*p6 + 2.174395E-14*p5 - 6.012288E-11*p4 + 8.772790E-08*p3 - 7.113554E-05*p2 + 3.029285E-02*p1 - 5.237677E+00;
138 }
139// a3
140 if (plab < 650.) {
141 a3=3.783071E-17*p6 - 1.151454E-13*p5 + 1.357165E-10*p4 - 8.036891E-08*p3 + 2.572396E-05*p2 - 4.245566E-03*p1 + 2.832772E-01;
142 }
143 else {
144 a3=-5.063316E-18*p6 + 3.223757E-14*p5 - 8.435635E-11*p4 + 1.159487E-07*p3 - 8.812510E-05*p2 + 3.500692E-02*p1 - 5.624556E+00;
145 }
146// a2
147 if (plab < 500.) {
148 a2=-6.085067E-14*p5 + 1.354078E-10*p4 - 1.124158E-07*p3 + 4.292106E-05*p2 - 7.218145E-03*p1 + 4.584962E-01;
149 }
150 else if (plab < 750.) {
151 a2= 9.512730E-11*p4 - 2.362724E-07*p3 + 2.171883E-04*p2 - 8.742722E-02*p1 + 1.309433E+01;
152 }
153 else {
154 a2=-4.228889E-18*p6 + 2.798222E-14*p5 - 7.640831E-11*p4 + 1.100124E-07*p3 - 8.778573E-05*p2 + 3.652772E-02*p1 - 6.025497E+00;
155 }
156// a1
157 if (plab < 500.) {
158 a1=-1.524408E-14*p5 + 3.007021E-11*p4 - 2.129570E-08*p3 + 5.607250E-06*p2 - 3.001598E-04*p1 + 8.701280E-04;
159 }
160 else if (plab < 750.) {
161 a1=-3.255396E-11*p4 + 8.168681E-08*p3 - 7.447474E-05*p2 + 2.917630E-02*p1 - 4.152037E+00;
162 }
163 else {
164 a1=9.964504E-19*p6 - 6.380168E-15*p5 + 1.638691E-11*p4 - 2.107063E-08*p3 + 1.347462E-05*p2 - 3.318304E-03*p1 - 5.030932E-02;
165 }
166// a0
167 a0=-3.220143E-17*p6 + 1.789654E-13*p5 - 3.912863E-10*p4 + 4.181510E-07*p3 - 2.147259E-04*p2 + 3.856266E-02*p1 + 2.609971E+00;
168
169 G4double interg1=2.*(a6/7. + a4/5. + a2/3. + a0); // (integral to normalize)
170 G4double f1=(a6+a5+a4+a3+a2+a1+a0)/interg1; // (Max normalized)
171
172 G4int passe1=0;
173 while (passe1==0) {
174 // Sample x from -1 to 1
175 x1=Random::shoot();
176 if (Random::shoot() > 0.5) x1=-x1;
177
178 // Sample u from 0 to 1
179 u1=Random::shoot();
180 fteta=(a6*x1*x1*x1*x1*x1*x1+a5*x1*x1*x1*x1*x1+a4*x1*x1*x1*x1+a3*x1*x1*x1+a2*x1*x1+a1*x1+a0)/interg1;
181 // The condition
182 if (u1*f1 < fteta) {
183 teta=std::acos(x1);
184 // std::cout << x1 << " " << fteta << " "<< f1/interg1 << " " << u1 << " " << interg1 << std::endl;
185 passe1=1;
186 }
187 }
188
189 fi=(2.0*pi)*Random::shoot();
190
191 ThreeVector mom_nucleon1(
192 pn*std::sin(teta)*std::cos(fi),
193 pn*std::sin(teta)*std::sin(fi),
194 pn*std::cos(teta)
195 );
196
197 mom_nucleon = -mom_nucleon1 ;
198
199 }
200
201 nucleon->setMomentum(mom_nucleon);
202 eta->setMomentum(-mom_nucleon);
203
204 fs->addModifiedParticle(nucleon);
205 fs->addModifiedParticle(eta);
206
207 }
208}
a0
const G4double a0
Definition: G4IonCoulombCrossSection.cc:68
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
G4INCL::EtaNElasticChannel::EtaNElasticChannel
EtaNElasticChannel(Particle *, Particle *)
Definition: G4INCLEtaNElasticChannel.cc:47
G4INCL::FinalState::addModifiedParticle
void addModifiedParticle(Particle *p)
Definition: G4INCLFinalState.cc:60
G4INCL::Particle::getEnergy
G4double getEnergy() const
Definition: G4INCLParticle.hh:774
G4INCL::Particle::setMomentum
virtual void setMomentum(const G4INCL::ThreeVector &momentum)
Definition: G4INCLParticle.hh:812
G4INCL::Particle::setEnergy
void setEnergy(G4double energy)
Definition: G4INCLParticle.hh:790
G4INCL::Particle::getMass
G4double getMass() const
Get the cached particle mass.
Definition: G4INCLParticle.hh:455
G4INCL::Particle::isNucleon
G4bool isNucleon() const
Definition: G4INCLParticle.hh:307
G4INCL::KinematicsUtils::momentumInLab
G4double momentumInLab(Particle const *const p1, Particle const *const p2)
gives the momentum in the lab frame of two particles.
Definition: G4INCLKinematicsUtils.cc:135
G4INCL::Random::normVector
ThreeVector normVector(G4double norm=1.)
Definition: G4INCLRandom.cc:162
G4INCL::Random::shoot
G4double shoot()
Definition: G4INCLRandom.cc:93