| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 H
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
45 namespace G4INCL {
46
47 EtaNElasticChannel::EtaNElasticChannel(Par
48 : particle1(p1), particle2(p2)
49 {
50
51 }
52
53 EtaNElasticChannel::~EtaNElasticChannel(){
54
55 }
56
57 void EtaNElasticChannel::fillFinalState(Fi
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::momen
69
70 G4double sh=nucleon->getEnergy()+eta
71 G4double mn=nucleon->getMass();
72 G4double me=eta->getMass();
73 G4double en=(sh*sh+mn*mn-me*me)/(2*s
74 nucleon->setEnergy(en);
75 G4double ee=std::sqrt(en*en-mn*mn+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.; // n
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-1
118 }
119 else {
120 a6=6.143615E-20*p6 - 3.157181E-1
121 }
122 // a5
123 if (plab < 650.) {
124 a5=-9.021076E-18*p6 + 2.176771E-
125 }
126 else if (plab < 950.){
127 a5=4.424756E-18*p6 - 1.756295E-1
128 }
129 else {
130 a5=2.209585E-19*p6 - 1.546647E-1
131 }
132 // a4
133 if (plab < 700.) {
134 a4=4.826684E-17*p6 - 1.534471E-1
135 }
136 else {
137 a4=-3.245143E-18*p6 + 2.174395E-
138 }
139 // a3
140 if (plab < 650.) {
141 a3=3.783071E-17*p6 - 1.151454E-1
142 }
143 else {
144 a3=-5.063316E-18*p6 + 3.223757E-
145 }
146 // a2
147 if (plab < 500.) {
148 a2=-6.085067E-14*p5 + 1.354078E-
149 }
150 else if (plab < 750.) {
151 a2= 9.512730E-11*p4 - 2.362724E-
152 }
153 else {
154 a2=-4.228889E-18*p6 + 2.798222E-
155 }
156 // a1
157 if (plab < 500.) {
158 a1=-1.524408E-14*p5 + 3.007021E-
159 }
160 else if (plab < 750.) {
161 a1=-3.255396E-11*p4 + 8.168681E-
162 }
163 else {
164 a1=9.964504E-19*p6 - 6.380168E-1
165 }
166 // a0
167 a0=-3.220143E-17*p6 + 1.789654E-
168
169 G4double interg1=2.*(a6/7. + a4/5.
170 G4double f1=(a6+a5+a4+a3+a2+a1+a0)
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=-x
177
178 // Sample u from 0 to 1
179 u1=Random::shoot();
180 fteta=(a6*x1*x1*x1*x1*x1*x1+a5*x
181 // The condition
182 if (u1*f1 < fteta) {
183 teta=std::acos(x1);
184 // std::cout << x1 <<
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 }
209