| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France 27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium 28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France 29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France 31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France
33 // 33 //
34 #define INCLXX_IN_GEANT4_MODE 1 34 #define INCLXX_IN_GEANT4_MODE 1
35 35
36 #include "globals.hh" 36 #include "globals.hh"
37 37
38 #include "G4INCLEtaNElasticChannel.hh" 38 #include "G4INCLEtaNElasticChannel.hh"
39 #include "G4INCLKinematicsUtils.hh" 39 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLBinaryCollisionAvatar.hh" 40 #include "G4INCLBinaryCollisionAvatar.hh"
41 #include "G4INCLRandom.hh" 41 #include "G4INCLRandom.hh"
42 #include "G4INCLGlobals.hh" 42 #include "G4INCLGlobals.hh"
43 #include "G4INCLLogger.hh" 43 #include "G4INCLLogger.hh"
44 44
45 namespace G4INCL { 45 namespace G4INCL {
46 46
47 EtaNElasticChannel::EtaNElasticChannel(Par 47 EtaNElasticChannel::EtaNElasticChannel(Particle *p1, Particle *p2)
48 : particle1(p1), particle2(p2) 48 : particle1(p1), particle2(p2)
49 { 49 {
50 50
51 } 51 }
52 52
53 EtaNElasticChannel::~EtaNElasticChannel(){ 53 EtaNElasticChannel::~EtaNElasticChannel(){
54 54
55 } 55 }
56 56
57 void EtaNElasticChannel::fillFinalState(Fi 57 void EtaNElasticChannel::fillFinalState(FinalState *fs) {
58 Particle * nucleon; 58 Particle * nucleon;
59 Particle * eta; 59 Particle * eta;
60 if(particle1->isNucleon()) { 60 if(particle1->isNucleon()) {
61 nucleon = particle1; 61 nucleon = particle1;
62 eta = particle2; 62 eta = particle2;
63 } else { 63 } else {
64 nucleon = particle2; 64 nucleon = particle2;
65 eta = particle1; 65 eta = particle1;
66 } 66 }
67 67
68 G4double plab=KinematicsUtils::momen 68 G4double plab=KinematicsUtils::momentumInLab(particle1, particle2);
69 69
70 G4double sh=nucleon->getEnergy()+eta 70 G4double sh=nucleon->getEnergy()+eta->getEnergy();
71 G4double mn=nucleon->getMass(); 71 G4double mn=nucleon->getMass();
72 G4double me=eta->getMass(); 72 G4double me=eta->getMass();
73 G4double en=(sh*sh+mn*mn-me*me)/(2*s 73 G4double en=(sh*sh+mn*mn-me*me)/(2*sh);
74 nucleon->setEnergy(en); 74 nucleon->setEnergy(en);
75 G4double ee=std::sqrt(en*en-mn*mn+me 75 G4double ee=std::sqrt(en*en-mn*mn+me*me);
76 eta->setEnergy(ee); 76 eta->setEnergy(ee);
77 G4double pn=std::sqrt(en*en-mn*mn); 77 G4double pn=std::sqrt(en*en-mn*mn);
78 78
79 ThreeVector mom_nucleon; 79 ThreeVector mom_nucleon;
80 80
81 if (plab < 250.) { 81 if (plab < 250.) {
82 // Isotropy 82 // Isotropy
83 mom_nucleon = Random::normVector(pn) 83 mom_nucleon = Random::normVector(pn);
84 } 84 }
85 85
86 // From Kamano 86 // From Kamano
87 else { 87 else {
88 88
89 const G4double pi=std::acos(-1.0); 89 const G4double pi=std::acos(-1.0);
90 G4double x1; 90 G4double x1;
91 G4double u1; 91 G4double u1;
92 G4double fteta; 92 G4double fteta;
93 G4double teta; 93 G4double teta;
94 G4double fi; 94 G4double fi;
95 95
96 G4double a0; 96 G4double a0;
97 G4double a1; 97 G4double a1;
98 G4double a2; 98 G4double a2;
99 G4double a3; 99 G4double a3;
100 G4double a4; 100 G4double a4;
101 G4double a5; 101 G4double a5;
102 G4double a6; 102 G4double a6;
103 103
104 if (plab > 1400.) plab=1400.; // n 104 if (plab > 1400.) plab=1400.; // no information on angular distributions above plab=1400 MeV
105 G4double p6=std::pow(plab, 6); 105 G4double p6=std::pow(plab, 6);
106 G4double p5=std::pow(plab, 5); 106 G4double p5=std::pow(plab, 5);
107 G4double p4=std::pow(plab, 4); 107 G4double p4=std::pow(plab, 4);
108 G4double p3=std::pow(plab, 3); 108 G4double p3=std::pow(plab, 3);
109 G4double p2=std::pow(plab, 2); 109 G4double p2=std::pow(plab, 2);
110 G4double p1=plab; 110 G4double p1=plab;
111 111
112 // a6 112 // a6
113 if (plab < 300.) { 113 if (plab < 300.) {
114 a6=-8.384000E-08*p1 - 1.15452E-04; 114 a6=-8.384000E-08*p1 - 1.15452E-04;
115 } 115 }
116 else if (plab < 500.){ 116 else if (plab < 500.){
117 a6=1.593966E-13*p4 - 2.619560E-1 117 a6=1.593966E-13*p4 - 2.619560E-10*p3 + 1.564701E-07*p2 - 3.986627E-05*p1 + 3.622575E-03;
118 } 118 }
119 else { 119 else {
120 a6=6.143615E-20*p6 - 3.157181E-1 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 } 121 }
122 // a5 122 // a5
123 if (plab < 650.) { 123 if (plab < 650.) {
124 a5=-9.021076E-18*p6 + 2.176771E- 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 } 125 }
126 else if (plab < 950.){ 126 else if (plab < 950.){
127 a5=4.424756E-18*p6 - 1.756295E-1 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 } 128 }
129 else { 129 else {
130 a5=2.209585E-19*p6 - 1.546647E-1 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 } 131 }
132 // a4 132 // a4
133 if (plab < 700.) { 133 if (plab < 700.) {
134 a4=4.826684E-17*p6 - 1.534471E-1 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 } 135 }
136 else { 136 else {
137 a4=-3.245143E-18*p6 + 2.174395E- 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 } 138 }
139 // a3 139 // a3
140 if (plab < 650.) { 140 if (plab < 650.) {
141 a3=3.783071E-17*p6 - 1.151454E-1 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 } 142 }
143 else { 143 else {
144 a3=-5.063316E-18*p6 + 3.223757E- 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 } 145 }
146 // a2 146 // a2
147 if (plab < 500.) { 147 if (plab < 500.) {
148 a2=-6.085067E-14*p5 + 1.354078E- 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 } 149 }
150 else if (plab < 750.) { 150 else if (plab < 750.) {
151 a2= 9.512730E-11*p4 - 2.362724E- 151 a2= 9.512730E-11*p4 - 2.362724E-07*p3 + 2.171883E-04*p2 - 8.742722E-02*p1 + 1.309433E+01;
152 } 152 }
153 else { 153 else {
154 a2=-4.228889E-18*p6 + 2.798222E- 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 } 155 }
156 // a1 156 // a1
157 if (plab < 500.) { 157 if (plab < 500.) {
158 a1=-1.524408E-14*p5 + 3.007021E- 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 } 159 }
160 else if (plab < 750.) { 160 else if (plab < 750.) {
161 a1=-3.255396E-11*p4 + 8.168681E- 161 a1=-3.255396E-11*p4 + 8.168681E-08*p3 - 7.447474E-05*p2 + 2.917630E-02*p1 - 4.152037E+00;
162 } 162 }
163 else { 163 else {
164 a1=9.964504E-19*p6 - 6.380168E-1 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 } 165 }
166 // a0 166 // a0
167 a0=-3.220143E-17*p6 + 1.789654E- 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 168
169 G4double interg1=2.*(a6/7. + a4/5. 169 G4double interg1=2.*(a6/7. + a4/5. + a2/3. + a0); // (integral to normalize)
170 G4double f1=(a6+a5+a4+a3+a2+a1+a0) 170 G4double f1=(a6+a5+a4+a3+a2+a1+a0)/interg1; // (Max normalized)
171 171
172 G4int passe1=0; 172 G4int passe1=0;
173 while (passe1==0) { 173 while (passe1==0) {
174 // Sample x from -1 to 1 174 // Sample x from -1 to 1
175 x1=Random::shoot(); 175 x1=Random::shoot();
176 if (Random::shoot() > 0.5) x1=-x 176 if (Random::shoot() > 0.5) x1=-x1;
177 177
178 // Sample u from 0 to 1 178 // Sample u from 0 to 1
179 u1=Random::shoot(); 179 u1=Random::shoot();
180 fteta=(a6*x1*x1*x1*x1*x1*x1+a5*x 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 181 // The condition
182 if (u1*f1 < fteta) { 182 if (u1*f1 < fteta) {
183 teta=std::acos(x1); 183 teta=std::acos(x1);
184 // std::cout << x1 << 184 // std::cout << x1 << " " << fteta << " "<< f1/interg1 << " " << u1 << " " << interg1 << std::endl;
185 passe1=1; 185 passe1=1;
186 } 186 }
187 } 187 }
188 188
189 fi=(2.0*pi)*Random::shoot(); 189 fi=(2.0*pi)*Random::shoot();
190 190
191 ThreeVector mom_nucleon1( 191 ThreeVector mom_nucleon1(
192 pn*std::sin(teta)*std::cos(fi), << 192 pn*std::sin(teta)*std::cos(fi),
193 pn*std::sin(teta)*std::sin(fi), << 193 pn*std::sin(teta)*std::sin(fi),
194 pn*std::cos(teta) << 194 pn*std::cos(teta)
195 ); << 195 );
196 196
197 mom_nucleon = -mom_nucleon1 ; 197 mom_nucleon = -mom_nucleon1 ;
198 198
199 } 199 }
200 200
201 nucleon->setMomentum(mom_nucleon); 201 nucleon->setMomentum(mom_nucleon);
202 eta->setMomentum(-mom_nucleon); 202 eta->setMomentum(-mom_nucleon);
203 203
204 fs->addModifiedParticle(nucleon); 204 fs->addModifiedParticle(nucleon);
205 fs->addModifiedParticle(eta); 205 fs->addModifiedParticle(eta);
206 206
207 } 207 }
208 } 208 }
209 209