| Geant4 Cross Reference |
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 "G4INCLNKbToNKbChannel.hh"
39 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLBinaryCollisionAvatar.hh"
41 #include "G4INCLRandom.hh"
42 #include "G4INCLGlobals.hh"
43 #include "G4INCLLogger.hh"
44 #include <algorithm>
45 #include "G4INCLPhaseSpaceGenerator.hh"
46
47 namespace G4INCL {
48
49 NKbToNKbChannel::NKbToNKbChannel(Particle *p1, Particle *p2)
50 : particle1(p1), particle2(p2)
51 {}
52
53 NKbToNKbChannel::~NKbToNKbChannel(){}
54
55 void NKbToNKbChannel::fillFinalState(FinalState *fs) {
56
57 Particle *nucleon;
58 Particle *kaon;
59
60 if(particle1->isNucleon()){
61 nucleon = particle1;
62 kaon = particle2;
63 }
64 else{
65 nucleon = particle2;
66 kaon = particle1;
67 }
68
69 // assert((ParticleTable::getIsospin(nucleon->getType()) + ParticleTable::getIsospin(kaon->getType())) == 0);
70
71 ThreeVector mom_kaon = KaonMomentum(kaon,nucleon);
72
73 if(kaon->getType() == KZeroBar){
74 nucleon->setType(Proton);
75 kaon->setType(KMinus);
76 }
77 else{
78 nucleon->setType(Neutron);
79 kaon->setType(KZeroBar);
80 }
81
82 G4double norm = KinematicsUtils::momentumInCM(kaon,nucleon);
83
84 kaon->setMomentum(mom_kaon*norm);
85 nucleon->setMomentum(-mom_kaon*norm);
86
87 kaon->adjustEnergyFromMomentum();
88 nucleon->adjustEnergyFromMomentum();
89
90
91 fs->addModifiedParticle(nucleon);
92 fs->addModifiedParticle(kaon);
93
94 }
95
96 ThreeVector NKbToNKbChannel::KaonMomentum(Particle const * const kaon, Particle const * const nucleon){
97
98 const G4double pLab = KinematicsUtils::momentumInLab(kaon,nucleon);
99
100 if(pLab < 235.) return Random::normVector(); // isotropic
101
102 G4double cos_theta = 1.;
103 G4double sin_theta = 0.;
104 const G4double cos_phi = std::cos(Random::shoot()*Math::twoPi);
105 const G4double sin_phi = std::sqrt(1-cos_phi*cos_phi);
106
107 const G4double x = kaon->getMomentum().getX();
108 const G4double y = kaon->getMomentum().getY();
109 const G4double z = kaon->getMomentum().getZ();
110
111 const G4double r = std::sqrt(x*x+y*y+z*z);
112 const G4double rho = std::sqrt(x*x+y*y);
113
114 if(pLab >= 1355.){
115 const G4double b = 12. * pLab/2375.; // correspond to the forward slope description at 2375 MeV/c in K- p elastic
116 cos_theta = std::log(Random::shoot()*(std::exp(b)-std::exp(-b))+std::exp(-b))/b;
117 sin_theta = std::sqrt(1-cos_theta*cos_theta);
118
119 }
120 else{
121 const G4double Legendre_coef[225][9] = {
122 {235,-0.30755,0,-0.04859,-0.01348,-9e-05,5e-05,0,0},
123 {240,-0.3041,0,-0.03549,-0.01096,-8e-05,5e-05,0,0},
124 {245,-0.30451,0,-0.02232,-0.00843,-8e-05,4e-05,0,0},
125 {250,-0.31028,0,-0.00899,-0.00585,-7e-05,4e-05,0,0},
126 {255,-0.32259,0,0.00462,-0.00319,-6e-05,4e-05,0,0},
127 {260,-0.34185,0,0.01859,-0.00041,-5e-05,3e-05,0,0},
128 {265,-0.36746,0,0.03305,0.00251,-4e-05,3e-05,0,0},
129 {270,-0.39766,0,0.04809,0.00563,-3e-05,2e-05,0,0},
130 {275,-0.42975,0,0.0638,0.00898,-2e-05,2e-05,0,0},
131 {280,-0.46069,1e-05,0.08023,0.0126,0,1e-05,0,0},
132 {285,-0.48736,5e-05,0.0973,0.0165,1e-05,0,0,0},
133 {290,-0.50729,2e-04,0.11487,0.02069,4e-05,-1e-05,0,0},
134 {295,-0.51855,0.00063,0.13282,0.02519,6e-05,-3e-05,0,0},
135 {300,-0.52024,0.00134,0.15096,0.02998,9e-05,-4e-05,0,0},
136 {305,-0.51207,0.00124,0.16892,0.03503,0.00013,-6e-05,0,0},
137 {310,-0.49468,-0.002,0.18618,0.04031,0.00016,-9e-05,0,0},
138 {315,-0.46892,-0.00933,0.20212,0.04568,0.00021,-0.00011,0,0},
139 {320,-0.43609,-0.01586,0.21604,0.05101,0.00026,-0.00014,0,0},
140 {325,-0.3973,-0.01303,0.22726,0.05617,0.00031,-0.00017,0,0},
141 {330,-0.35379,0.00167,0.23511,0.06107,0.00037,-0.00021,0,0},
142 {335,-0.3067,0.0219,0.23911,0.06576,0.00044,-0.00024,0,0},
143 {340,-0.25739,0.04678,0.23878,0.07033,5e-04,-0.00028,0,0},
144 {345,-0.20755,0.09094,0.23367,0.07488,0.00057,-0.00032,0,0},
145 {350,-0.15908,0.17482,0.22334,0.0795,0.00064,-0.00037,0,0},
146 {355,-0.11446,0.30682,0.20753,0.08421,0.00071,-0.00041,0,0},
147 {360,-0.07625,0.48075,0.1861,0.08906,0.00078,-0.00045,0,0},
148 {365,-0.04747,0.68345,0.1593,0.09421,0.00083,-0.00049,0,0},
149 {370,-0.03064,0.90388,0.1275,0.09983,0.00088,-0.00052,0,0},
150 {375,-0.02812,1.13397,0.09104,0.10608,0.00091,-0.00054,0,0},
151 {380,-0.04087,1.3636,0.05039,0.11311,0.00092,-0.00056,0,0},
152 {385,-0.06879,1.57657,0.00644,0.12093,9e-04,-0.00056,0,0},
153 {390,-0.10972,1.75346,-0.03978,0.1295,0.00085,-0.00054,0,0},
154 {395,-0.16022,1.87981,-0.08707,0.13859,0.00076,-5e-04,0,0},
155 {400,-0.21553,1.9509,-0.1342,0.14795,0.00063,-0.00043,0,0},
156 {405,-0.27043,1.97287,-0.17995,0.15732,0.00043,-0.00034,0,0},
157 {410,-0.32009,1.95458,-0.22323,0.16643,0.00017,-2e-04,0,0},
158 {415,-0.36095,1.89945,-0.26354,0.1749,-0.00016,-2e-05,0,0},
159 {420,-0.3912,1.80489,-0.30052,0.18232,-0.00058,2e-04,0,0},
160 {425,-0.4106,1.68228,-0.3338,0.18828,-0.0011,0.00048,0,0},
161 {430,-0.42021,1.56358,-0.36315,0.19246,-0.00172,0.00082,0,0},
162 {435,-0.42128,1.47378,-0.38888,0.19478,-0.00244,0.00123,0,0},
163 {440,-0.415,1.42688,-0.41138,0.19524,-0.00329,0.0017,0,0},
164 {445,-0.40201,1.39998,-0.43103,0.19376,-0.00425,0.00225,0,0},
165 {450,-0.38247,1.36801,-0.44826,0.19036,-0.00534,0.00286,0,0},
166 {455,-0.35589,1.33185,-0.46368,0.18542,-0.00653,0.00355,0,0},
167 {460,-0.3216,1.29215,-0.47787,0.17941,-0.00782,0.0043,0,0},
168 {465,-0.27917,1.2499,-0.49102,0.1727,-0.00918,0.0051,0,0},
169 {470,-0.22973,1.20597,-0.50321,0.16567,-0.0106,0.00594,0,0},
170 {475,-0.17713,1.16167,-0.5145,0.15865,-0.01204,0.00682,0,0},
171 {480,-0.12841,1.11813,-0.52493,0.15196,-0.01349,0.00771,0,0},
172 {485,-0.09149,1.07666,-0.53436,0.14575,-0.01489,0.0086,0,0},
173 {490,-0.07028,1.03839,-0.5426,0.14013,-0.01618,0.00944,0,0},
174 {495,-0.06296,1.00433,-0.54947,0.13518,-0.0173,0.01021,0,0},
175 {500,-0.06366,0.97503,-0.55477,0.13101,-0.01819,0.01087,0,0},
176 {505,-0.06664,0.95085,-0.55833,0.12773,-0.0188,0.0114,0,0},
177 {510,-0.06771,0.93169,-0.55996,0.12544,-0.01906,0.01175,0,0},
178 {515,-0.06476,0.91724,-0.55947,0.12423,-0.01891,0.01191,0,0},
179 {520,-0.05705,0.90697,-0.55703,0.12437,-0.01831,0.01186,0,0},
180 {525,-0.04488,0.90021,-0.55318,0.12594,-0.01736,0.01179,0,0},
181 {530,-0.02927,0.8963,-0.54843,0.12893,-0.01624,0.01196,0,0},
182 {535,-0.01183,0.8945,-0.5433,0.13334,-0.01512,0.0126,0,0},
183 {540,0.00596,0.89424,-0.53831,0.13906,-0.01414,0.01394,0,0},
184 {545,0.02306,0.89492,-0.5339,0.1458,-0.0134,0.01611,0,0},
185 {550,0.039,0.89617,-0.53048,0.15329,-0.01298,0.01915,0,0},
186 {555,0.05391,0.89758,-0.52847,0.16125,-0.01294,0.02316,0,0},
187 {560,0.06806,0.89897,-0.52824,0.16947,-0.01334,0.02816,0,0},
188 {565,0.08166,0.90012,-0.52974,0.17794,-0.01416,0.03398,0,0},
189 {570,0.09478,0.90105,-0.53283,0.18651,-0.01529,0.0403,0,0},
190 {575,0.10746,0.9017,-0.53736,0.19494,-0.01659,0.04682,0,0},
191 {580,0.11985,0.90216,-0.54318,0.20297,-0.01795,0.05322,0,0},
192 {585,0.13218,0.90248,-0.55002,0.21039,-0.01909,0.05908,0,0},
193 {590,0.14481,0.90282,-0.55752,0.21698,-0.01964,0.06395,0,0},
194 {595,0.15769,0.90328,-0.56531,0.22255,-0.01921,0.06734,0,0},
195 {600,0.17046,0.90397,-0.57302,0.2269,-0.01743,0.0688,0,0},
196 {605,0.18205,0.90496,-0.58026,0.23009,-0.01427,0.06815,0,0},
197 {610,0.19103,0.90635,-0.58673,0.23229,-0.00999,0.06549,0,0},
198 {615,0.19589,0.90814,-0.59209,0.2337,-0.00486,0.06093,0,0},
199 {620,0.19545,0.91035,-0.59603,0.2345,0.00083,0.05457,0,0},
200 {625,0.18873,0.91297,-0.59864,0.23485,0.00685,0.04682,0,0},
201 {630,0.17514,0.91598,-0.60046,0.2349,0.01298,0.03844,0,0},
202 {635,0.15461,0.91935,-0.60188,0.23486,0.0191,0.03017,0,0},
203 {640,0.12791,0.92305,-0.60283,0.23506,0.0253,0.02263,0,0},
204 {645,0.09694,0.92707,-0.60322,0.23576,0.03151,0.01633,0,0},
205 {650,0.06435,0.9314,-0.60304,0.23717,0.03758,0.01172,0,0},
206 {655,0.03278,0.93603,-0.60224,0.23951,0.04331,0.00924,0,0},
207 {660,0.00403,0.94096,-0.60081,0.24298,0.04854,0.00932,0,0},
208 {665,-0.02138,0.9462,-0.59882,0.24774,0.05304,0.01201,0,0},
209 {670,-0.04399,0.95174,-0.59638,0.25384,0.05655,0.01685,0,0},
210 {675,-0.06489,0.95757,-0.59359,0.26128,0.05879,0.02333,0,0},
211 {680,-0.08497,0.96367,-0.59053,0.27011,0.05955,0.03096,-1e-05,0},
212 {685,-0.10471,0.97,-0.58712,0.2802,0.0589,0.03927,-1e-05,0},
213 {690,-0.12371,0.97649,-0.58315,0.29136,0.05721,0.0478,-1e-05,0},
214 {695,-0.14097,0.98305,-0.57828,0.3031,0.05487,0.05603,-1e-05,0},
215 {700,-0.15469,0.98957,-0.57209,0.31472,0.05231,0.06337,-1e-05,0},
216 {705,-0.16271,0.99591,-0.56443,0.3256,0.04979,0.06932,-1e-05,0},
217 {710,-0.16255,1.00194,-0.55574,0.33535,0.04732,0.07359,-2e-05,0},
218 {715,-0.15199,1.00752,-0.54657,0.34363,0.04488,0.07591,-2e-05,0},
219 {720,-0.12951,1.01251,-0.53743,0.3499,0.04245,0.07587,-2e-05,0},
220 {725,-0.09529,1.0168,-0.52879,0.35376,0.03995,0.07327,-2e-05,0},
221 {730,-0.05159,1.0203,-0.5212,0.35562,0.03722,0.06853,-2e-05,0},
222 {735,-0.00245,1.02293,-0.51519,0.35595,0.0341,0.06213,-2e-05,0},
223 {740,0.04772,1.02468,-0.51129,0.35522,0.03044,0.05455,-2e-05,0},
224 {745,0.09584,1.02553,-0.51005,0.35387,0.02597,0.04627,-2e-05,0},
225 {750,0.14086,1.02552,-0.51202,0.35229,0.02027,0.03776,-1e-05,0},
226 {755,0.18338,1.02466,-0.51778,0.35084,0.01302,0.02948,-1e-05,0},
227 {760,0.22444,1.023,-0.52787,0.34965,0.00409,0.02179,0,0},
228 {765,0.2645,1.0206,-0.54268,0.34874,-0.00656,0.01494,1e-05,0},
229 {770,0.30276,1.01748,-0.56211,0.34808,-0.01882,0.00896,2e-05,0},
230 {775,0.33729,1.01369,-0.58582,0.34762,-0.03243,0.00376,4e-05,0},
231 {780,0.36532,1.00926,-0.61343,0.34741,-0.04713,-8e-04,6e-05,0},
232 {785,0.38445,1.00421,-0.64453,0.34758,-0.06254,-0.00496,8e-05,0},
233 {790,0.39323,0.99857,-0.67869,0.34828,-0.07831,-0.00897,0.00011,0},
234 {795,0.39129,0.99243,-0.71549,0.34967,-0.09409,-0.01306,0.00014,0},
235 {800,0.37889,0.98587,-0.75456,0.35269,-0.10967,-0.01765,0.00017,0},
236 {805,0.35624,0.97903,-0.79561,0.35848,-0.12507,-0.023,0.00021,0},
237 {810,0.32292,0.97211,-0.83834,0.36717,-0.14022,-0.02894,0.00025,0},
238 {815,0.27791,0.96533,-0.88246,0.37885,-0.15509,-0.03525,0.00029,0},
239 {820,0.22029,0.95894,-0.92768,0.39359,-0.16961,-0.04173,0.00034,0},
240 {825,0.15048,0.9532,-0.9737,0.41145,-0.18373,-0.04817,0.00039,0},
241 {830,0.07108,0.94826,-1.02023,0.43252,-0.1974,-0.05435,0.00044,0},
242 {835,-0.0136,0.94426,-1.06697,0.45687,-0.21057,-0.06007,0.00049,0},
243 {840,-0.09855,0.94118,-1.11361,0.48455,-0.22318,-0.06513,0.00054,0},
244 {845,-0.17951,0.9389,-1.15976,0.51525,-0.23515,-0.06938,0.00059,0},
245 {850,-0.25366,0.93727,-1.20484,0.54811,-0.24638,-0.07279,0.00064,0},
246 {855,-0.31965,0.93606,-1.24827,0.58227,-0.25672,-0.07535,0.00069,0},
247 {860,-0.37707,0.93506,-1.28948,0.61685,-0.26608,-0.07701,0.00072,0},
248 {865,-0.42612,0.93404,-1.3279,0.651,-0.27433,-0.07778,0.00075,0},
249 {870,-0.46711,0.93281,-1.36298,0.68418,-0.28136,-0.0779,0.00077,0},
250 {875,-0.50042,0.93129,-1.39424,0.7161,-0.28708,-0.07783,0.00078,0},
251 {880,-0.52649,0.92936,-1.42118,0.74648,-0.29139,-0.07803,0.00076,0},
252 {885,-0.54572,0.92699,-1.44339,0.77492,-0.29417,-0.07883,0.00072,0},
253 {890,-0.55824,0.92417,-1.46054,0.801,-0.29529,-0.08047,0.00065,0},
254 {895,-0.56399,0.92098,-1.47259,0.82492,-0.29497,-0.08289,0.00054,0},
255 {900,-0.56271,0.91754,-1.47975,0.84731,-0.29367,-0.08579,0.00039,0},
256 {905,-0.55466,0.91396,-1.48218,0.86857,-0.29166,-0.08892,0.00018,0},
257 {910,-0.54086,0.91044,-1.48007,0.88907,-0.28921,-0.09206,-8e-05,0},
258 {915,-0.52332,0.90717,-1.47358,0.90921,-0.28662,-0.09496,-0.00042,0},
259 {920,-0.5045,0.90438,-1.46288,0.92935,-0.28415,-0.09741,-0.00084,0},
260 {925,-0.48654,0.9023,-1.4482,0.94975,-0.28214,-0.09924,-0.00135,0},
261 {930,-0.47053,0.90114,-1.42984,0.97058,-0.281,-0.10034,-0.00195,0},
262 {935,-0.45659,0.90112,-1.40806,0.99201,-0.2811,-0.10059,-0.00266,0},
263 {940,-0.44423,0.90246,-1.38314,1.01421,-0.28285,-0.09988,-0.00348,0},
264 {945,-0.43287,0.90525,-1.35536,1.0371,-0.28653,-0.09816,-0.00441,0},
265 {950,-0.42193,0.90965,-1.32498,1.06041,-0.29239,-0.09538,-0.00544,0},
266 {955,-0.41093,0.91577,-1.29226,1.08389,-0.30066,-0.09153,-0.00659,0},
267 {960,-0.39924,0.92369,-1.25746,1.10725,-0.31158,-0.08656,-0.00783,0},
268 {965,-0.38604,0.93346,-1.22095,1.13008,-0.32524,-0.08039,-0.00918,0},
269 {970,-0.37027,0.94514,-1.18303,1.15197,-0.34144,-0.07311,-0.01058,0},
270 {975,-0.35077,0.95876,-1.14385,1.17271,-0.35991,-0.06496,-0.01199,0},
271 {980,-0.32717,0.97435,-1.10358,1.19206,-0.38038,-0.05616,-0.01336,0},
272 {985,-0.3004,0.99188,-1.0624,1.20979,-0.40261,-0.04696,-0.01464,0},
273 {990,-0.27312,1.01135,-1.02046,1.2257,-0.42632,-0.03758,-0.01578,0},
274 {995,-0.24834,1.03272,-0.97795,1.23954,-0.45126,-0.02826,-0.01674,0},
275 {1000,-0.22769,1.0559,-0.93501,1.2511,-0.47716,-0.01922,-0.01745,0},
276 {1005,-0.21032,1.0808,-0.89184,1.26017,-0.50377,-0.01069,-0.01787,0},
277 {1010,-0.19369,1.10728,-0.84866,1.26669,-0.53095,-0.00261,-0.0179,0},
278 {1015,-0.17436,1.13513,-0.80576,1.27066,-0.55859,0.00512,-0.01742,0},
279 {1020,-0.1494,1.16411,-0.7634,1.27208,-0.58657,0.01264,-0.01633,0},
280 {1025,-0.11749,1.19391,-0.72184,1.27097,-0.61479,0.02007,-0.01451,0},
281 {1030,-0.081,1.22417,-0.68137,1.26732,-0.64314,0.02755,-0.01185,0},
282 {1035,-0.04533,1.25445,-0.64225,1.26114,-0.6715,0.0352,-0.00824,0},
283 {1040,-0.01646,1.28424,-0.60476,1.25244,-0.69977,0.04316,-0.00356,0},
284 {1045,0.00361,1.31304,-0.56915,1.24121,-0.72785,0.05154,0.00231,0},
285 {1050,0.01552,1.34024,-0.53567,1.22751,-0.75561,0.06047,0.0094,0},
286 {1055,0.02103,1.36529,-0.50435,1.21159,-0.78295,0.07001,0.01752,0},
287 {1060,0.02215,1.38771,-0.47519,1.19376,-0.80975,0.08022,0.02639,0},
288 {1065,0.02047,1.40712,-0.4482,1.1743,-0.83589,0.09116,0.03575,0},
289 {1070,0.01735,1.4233,-0.42337,1.15353,-0.86126,0.10289,0.04532,0},
290 {1075,0.0137,1.43585,-0.40071,1.13173,-0.88574,0.11546,0.05485,0},
291 {1080,0.00983,1.44489,-0.38022,1.10922,-0.90921,0.12892,0.06405,0},
292 {1085,0.00516,1.45066,-0.3619,1.08628,-0.93157,0.14335,0.07265,0},
293 {1090,-0.00091,1.45292,-0.34574,1.0632,-0.95272,0.15877,0.08044,0},
294 {1095,-0.00921,1.4522,-0.33171,1.04013,-0.97273,0.17517,0.08734,0},
295 {1100,-0.01996,1.44884,-0.31977,1.01722,-0.99171,0.19252,0.09332,-2e-05},
296 {1105,-0.03239,1.44326,-0.30987,0.99459,-1.00975,0.2108,0.09835,-4e-05},
297 {1110,-0.0447,1.43582,-0.30197,0.97238,-1.02695,0.22997,0.1024,-6e-05},
298 {1115,-0.05508,1.42673,-0.29604,0.95074,-1.0434,0.25,0.10545,-9e-05},
299 {1120,-0.06264,1.41667,-0.29202,0.92979,-1.05922,0.27085,0.10746,-0.00015},
300 {1125,-0.06741,1.40585,-0.28988,0.90967,-1.0745,0.29251,0.10841,-0.00023},
301 {1130,-0.06996,1.39476,-0.28958,0.89052,-1.08934,0.31494,0.10826,-0.00036},
302 {1135,-0.07105,1.38376,-0.29107,0.87247,-1.10383,0.3381,0.107,-0.00054},
303 {1140,-0.0711,1.3732,-0.29431,0.85567,-1.11808,0.36197,0.10458,-0.00087},
304 {1145,-0.07085,1.36335,-0.29927,0.84024,-1.13219,0.38652,0.10099,-0.00128},
305 {1150,-0.07116,1.35452,-0.30589,0.82632,-1.14626,0.41172,0.09619,-0.00185},
306 {1155,-0.07304,1.34687,-0.31414,0.81405,-1.16038,0.43753,0.09015,-0.00263},
307 {1160,-0.07771,1.34061,-0.32397,0.80356,-1.17466,0.46393,0.08284,-0.00367},
308 {1165,-0.08664,1.3359,-0.33535,0.79499,-1.18919,0.49089,0.07424,-0.00506},
309 {1170,-0.10173,1.33284,-0.34823,0.78841,-1.20404,0.51834,0.06434,-0.00686},
310 {1175,-0.12509,1.33145,-0.36258,0.78356,-1.21914,0.54614,0.05318,-0.00925},
311 {1180,-0.15682,1.33177,-0.37838,0.78014,-1.23438,0.57409,0.04086,-0.0122},
312 {1185,-0.19464,1.33374,-0.39559,0.77782,-1.24965,0.60201,0.02744,-0.01586},
313 {1190,-0.23352,1.33727,-0.41419,0.7763,-1.26484,0.62972,0.01301,-0.02035},
314 {1195,-0.26876,1.34221,-0.43415,0.77525,-1.27985,0.65702,-0.00237,-0.02577},
315 {1200,-0.299,1.34841,-0.45544,0.77437,-1.29457,0.68375,-0.01862,-0.0322},
316 {1205,-0.32631,1.35562,-0.47804,0.77332,-1.30888,0.7097,-0.03565,-0.03976},
317 {1210,-0.35505,1.36358,-0.50188,0.77186,-1.32269,0.73474,-0.05339,-0.04851},
318 {1215,-0.38961,1.372,-0.52686,0.7699,-1.33587,0.75888,-0.07178,-0.0585},
319 {1220,-0.43176,1.38055,-0.55282,0.76744,-1.3483,0.78219,-0.09072,-0.06977},
320 {1225,-0.47936,1.38889,-0.57962,0.76446,-1.35986,0.80471,-0.11016,-0.08229},
321 {1230,-0.52635,1.39666,-0.60711,0.76094,-1.37043,0.82652,-0.13,-0.09604},
322 {1235,-0.56649,1.40351,-0.63516,0.75687,-1.37988,0.84766,-0.15019,-0.11094},
323 {1240,-0.59709,1.40914,-0.66361,0.75222,-1.3881,0.86819,-0.17064,-0.12681},
324 {1245,-0.61879,1.41323,-0.69232,0.74699,-1.39497,0.88817,-0.19128,-0.14349},
325 {1250,-0.63452,1.41555,-0.72116,0.74115,-1.40036,0.90767,-0.21203,-0.16074},
326 {1255,-0.6472,1.4159,-0.74996,0.73469,-1.40415,0.92673,-0.23283,-0.17827},
327 {1260,-0.65841,1.4142,-0.77859,0.7276,-1.40622,0.94542,-0.25358,-0.19575},
328 {1265,-0.66819,1.41043,-0.80691,0.71985,-1.40646,0.96379,-0.27423,-0.21283},
329 {1270,-0.67506,1.40465,-0.83476,0.71143,-1.40474,0.9819,-0.29469,-0.22916},
330 {1275,-0.67695,1.39703,-0.86202,0.70233,-1.40093,0.99981,-0.3149,-0.24439},
331 {1280,-0.67215,1.38781,-0.88852,0.69252,-1.39492,1.01758,-0.33477,-0.25814},
332 {1285,-0.65883,1.37728,-0.91412,0.68199,-1.38659,1.03527,-0.35423,-0.27016},
333 {1290,-0.63639,1.36579,-0.93873,0.67076,-1.37587,1.05293,-0.37323,-0.28021},
334 {1295,-0.60463,1.35371,-0.96237,0.65893,-1.3629,1.07059,-0.3918,-0.28812},
335 {1300,-0.56527,1.34141,-0.98514,0.64665,-1.34787,1.08829,-0.41,-0.29382},
336 {1305,-0.52229,1.32923,-1.0071,0.63407,-1.33099,1.10606,-0.4279,-0.29727},
337 {1310,-0.4806,1.31751,-1.02834,0.62132,-1.31244,1.12394,-0.44553,-0.29853},
338 {1315,-0.44384,1.3065,-1.04893,0.60855,-1.29241,1.14195,-0.46298,-0.29776},
339 {1320,-0.41298,1.29641,-1.06896,0.5959,-1.27111,1.16013,-0.48028,-0.29505},
340 {1325,-0.38708,1.28739,-1.0885,0.58352,-1.24873,1.17852,-0.4975,-0.29066},
341 {1330,-0.36429,1.27952,-1.10764,0.57154,-1.22546,1.19714,-0.5147,-0.28485},
342 {1335,-0.34347,1.27284,-1.12645,0.56012,-1.2015,1.21602,-0.53193,-0.27787},
343 {1340,-0.32464,1.26732,-1.14501,0.54939,-1.17704,1.23521,-0.54926,-0.26998},
344 {1345,-0.30854,1.2629,-1.1634,0.5395,-1.15228,1.25472,-0.56673,-0.26145},
345 {1350,-0.29584,1.25949,-1.18169,0.53059,-1.1274,1.27461,-0.58441,-0.2525},
346 {1355,-0.28655,1.25697,-1.19998,0.52281,-1.10262,1.29488,-0.60235,-0.2434}};
347
348 const G4int coef_ener = G4int((pLab-Legendre_coef[0][0])/5);
349 const G4double sup_ener = pLab/5. - coef_ener -Legendre_coef[0][0]/5;
350
351 // assert(pLab >= Legendre_coef[coef_ener][0] && pLab < Legendre_coef[coef_ener+1][0]);
352
353 // Legendre coefficient normalized
354 const G4double A0 = 1.;
355 const G4double A1 = (1-sup_ener)*Legendre_coef[coef_ener][1] + sup_ener*Legendre_coef[coef_ener+1][1];
356 const G4double A2 = (1-sup_ener)*Legendre_coef[coef_ener][2] + sup_ener*Legendre_coef[coef_ener+1][2];
357 const G4double A3 = (1-sup_ener)*Legendre_coef[coef_ener][3] + sup_ener*Legendre_coef[coef_ener+1][3];
358 const G4double A4 = (1-sup_ener)*Legendre_coef[coef_ener][4] + sup_ener*Legendre_coef[coef_ener+1][4];
359 const G4double A5 = (1-sup_ener)*Legendre_coef[coef_ener][5] + sup_ener*Legendre_coef[coef_ener+1][5];
360 const G4double A6 = (1-sup_ener)*Legendre_coef[coef_ener][6] + sup_ener*Legendre_coef[coef_ener+1][6];
361 const G4double A7 = (1-sup_ener)*Legendre_coef[coef_ener][7] + sup_ener*Legendre_coef[coef_ener+1][7];
362 const G4double A8 = (1-sup_ener)*Legendre_coef[coef_ener][8] + sup_ener*Legendre_coef[coef_ener+1][8];
363
364 // Theoritical max if all Ai > 0 (often the case)
365 const G4double A = std::fabs(A0) + std::fabs(A1) + std::fabs(A2) + std::fabs(A3) + std::fabs(A4) + std::fabs(A5) + std::fabs(A6) + std::fabs(A7) + std::fabs(A8);
366
367 G4bool success = false;
368 G4int maxloop = 0;
369
370 while(!success && maxloop < 1000){
371
372 cos_theta = Random::shoot()*2-1.; // not optimized
373
374 // Legendre Polynomial
375 G4double P0 = A0;
376 G4double P1 = A1*cos_theta;
377 G4double P2 = A2/2.*(3*std::pow(cos_theta,2)-1);
378 G4double P3 = A3/2.*(5*std::pow(cos_theta,3)-3*cos_theta);
379 G4double P4 = A4/8.*(35*std::pow(cos_theta,4)-30*std::pow(cos_theta,2)+3);
380 G4double P5 = A5/8.*(63*std::pow(cos_theta,5)-70*std::pow(cos_theta,3)+15*cos_theta);
381 G4double P6 = A6/16.*(231*std::pow(cos_theta,6)-315*std::pow(cos_theta,4)+105*std::pow(cos_theta,2)-5);
382 G4double P7 = A7/16.*(429*std::pow(cos_theta,7)-693*std::pow(cos_theta,5)+315*std::pow(cos_theta,3)-35*cos_theta);
383 G4double P8 = A8/128.*(6435*std::pow(cos_theta,8)-12012*std::pow(cos_theta,6)+6930*std::pow(cos_theta,4)-1260*std::pow(cos_theta,2)+35);
384
385 G4double P = (P0 + P1 + P2 + P3 + P4 + P5 + P6 + P7 + P8)/2; // /2 for the normalisation
386
387 if(Random::shoot()*A < P) success = true;
388 maxloop +=1 ;
389 if(maxloop==1000) cos_theta = std::log(Random::shoot()*(std::exp(10.)-std::exp(-10.))+std::exp(-10.))/10.; // if no success in 1E4 shoot, probably angulard distribution piked very forward
390 }
391 sin_theta = std::sqrt(1-cos_theta*cos_theta);
392 }
393
394 if(rho == 0) return ThreeVector(sin_theta*cos_phi,sin_theta*sin_phi,cos_theta);
395 // Rotation in the direction of the incident kaon
396 const G4double px = x/r*cos_theta - y/rho*sin_theta*cos_phi + z/r*x/rho*sin_theta*sin_phi;
397 const G4double py = y/r*cos_theta + x/rho*sin_theta*cos_phi + z/r*y/rho*sin_theta*sin_phi;
398 const G4double pz = z/r*cos_theta - rho/r*sin_theta*sin_phi;
399
400 return ThreeVector(px,py,pz);
401 }
402 }
403