ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLNpiToLKChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLNpiToLKChannel.cc
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 "G4INCLNpiToLKChannel.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  NpiToLKChannel::NpiToLKChannel(Particle *p1, Particle *p2)
50  : particle1(p1), particle2(p2)
51  {}
52 
53  NpiToLKChannel::~NpiToLKChannel(){}
54 
55  void NpiToLKChannel::fillFinalState(FinalState *fs) {
56 
57  Particle *nucleon;
58  Particle *pion;
59 
60  if(particle1->isNucleon()){
61  nucleon = particle1;
62  pion = particle2;
63  }
64  else{
65  nucleon = particle2;
66  pion = particle1;
67  }
68 
69  const G4int iso = ParticleTable::getIsospin(nucleon->getType()) + ParticleTable::getIsospin(pion->getType());
70 
71  ParticleType KaonType;
72 
73  if(iso == 1) KaonType = KPlus;
74  else if(iso == -1) KaonType = KZero;
75  else{
76  INCL_ERROR("NpiToLKChannel called with an inconsitant pair\n");
77  return;
78  }
79 
80  ThreeVector mom_kaon = KaonMomentum(pion,nucleon);
81 
82  nucleon->setType(Lambda);
83  pion->setType(KaonType);
84 
85  G4double norm = KinematicsUtils::momentumInCM(nucleon,pion);
86 
87  pion->setMomentum(mom_kaon*norm);
88  nucleon->setMomentum(-mom_kaon*norm);
89 
90  nucleon->adjustEnergyFromMomentum();
91  pion->adjustEnergyFromMomentum();
92 
93  //INCL_DEBUG("NpiToLK " << (pion->getMomentum().theta()) * 180. / G4INCL::Math::pi << '\n');
94 
95  fs->addModifiedParticle(nucleon);
96  fs->addModifiedParticle(pion);
97 
98  }
99 
100  ThreeVector NpiToLKChannel::KaonMomentum(Particle const * const pion, Particle const * const nucleon){
101 
102  const G4double pLab = KinematicsUtils::momentumInLab(pion,nucleon);
103 
104  if(pLab < 930.) return Random::normVector(); // isotropic
105 
106  G4double cos_theta = 1.;
107  G4double sin_theta = 0.;
108  const G4double cos_phi = std::cos(Random::shoot()*Math::twoPi);
109  const G4double sin_phi = std::sqrt(1-cos_phi*cos_phi);
110 
111  const G4double x = pion->getMomentum().getX();
112  const G4double y = pion->getMomentum().getY();
113  const G4double z = pion->getMomentum().getZ();
114 
115  const G4double r = std::sqrt(x*x+y*y+z*z);
116  const G4double rho = std::sqrt(x*x+y*y);
117 
118  if(pLab >= 2375.){
119  const G4double b = 12. * pLab/2375.; // correspond to the forward slope description at 2375 MeV/c
120  cos_theta = std::log(Random::shoot()*(std::exp(b)-std::exp(-b))+std::exp(-b))/b;
121  sin_theta = std::sqrt(1-cos_theta*cos_theta);
122 
123  }
124  else{
125  const G4double Legendre_coef[290][8] = {
126  {930,0.18739,-0.06,-0.00621,-2e-05,0,0,0},
127  {935,0.22915,-0.05104,-0.00013,-2e-05,0,0,0},
128  {940,0.27087,-0.04199,0.006,-2e-05,0,0,0},
129  {945,0.31247,-0.0327,0.01223,-2e-05,0,0,0},
130  {950,0.35388,-0.02306,0.01863,-2e-05,0,0,0},
131  {955,0.39503,-0.01292,0.02525,-2e-05,0,0,0},
132  {960,0.43586,-0.00217,0.03214,-3e-05,0,0,0},
133  {965,0.4763,0.00933,0.03938,-3e-05,0,0,0},
134  {970,0.51629,0.02171,0.04702,-3e-05,0,0,0},
135  {975,0.55575,0.03509,0.05511,-3e-05,0,0,0},
136  {980,0.59462,0.04962,0.06371,-3e-05,0,0,0},
137  {985,0.63291,0.06537,0.073,-3e-05,0,0,0},
138  {990,0.6706,0.08222,0.0833,-3e-05,0,0,0},
139  {995,0.70756,0.09998,0.09494,-2e-05,0,0,0},
140  {1000,0.74362,0.11842,0.10816,-2e-05,0,0,0},
141  {1005,0.77857,0.13733,0.12261,-2e-05,0,0,0},
142  {1010,0.8121,0.15648,0.13737,-1e-05,0,0,0},
143  {1015,0.84402,0.17556,0.15102,-1e-05,0,0,0},
144  {1020,0.87447,0.19424,0.16201,0,0,0,0},
145  {1025,0.90358,0.21231,0.16943,1e-05,0,0,0},
146  {1030,0.93118,0.22961,0.17264,2e-05,0,0,0},
147  {1035,0.95642,0.24598,0.17205,4e-05,0,0,0},
148  {1040,0.9784,0.26128,0.16958,6e-05,0,0,0},
149  {1045,0.99681,0.2754,0.16737,8e-05,0,0,0},
150  {1050,1.01213,0.2883,0.16575,1e-04,0,0,0},
151  {1055,1.02473,0.29992,0.16397,0.00012,0,0,0},
152  {1060,1.03467,0.31018,0.16082,0.00015,0,0,0},
153  {1065,1.0417,0.31907,0.15599,0.00018,0,0,0},
154  {1070,1.04567,0.32661,0.14947,0.00021,0,0,0},
155  {1075,1.04705,0.3329,0.14185,0.00024,0,0,0},
156  {1080,1.04645,0.33806,0.13383,0.00027,0,0,0},
157  {1085,1.04451,0.34222,0.12613,3e-04,0,0,0},
158  {1090,1.04182,0.34551,0.11946,0.00032,0,0,0},
159  {1095,1.03903,0.34804,0.11454,0.00034,0,0,0},
160  {1100,1.03662,0.34992,0.11184,0.00036,0,0,0},
161  {1105,1.03464,0.35117,0.11087,0.00037,0,0,0},
162  {1110,1.03303,0.3518,0.11087,0.00037,0,0,0},
163  {1115,1.0317,0.35178,0.1111,0.00037,0,0,0},
164  {1120,1.03059,0.35112,0.11082,0.00035,0,0,0},
165  {1125,1.02963,0.34981,0.10927,0.00031,0,0,0},
166  {1130,1.02873,0.34786,0.10573,0.00027,0,0,0},
167  {1135,1.02779,0.34528,0.09967,2e-04,0,0,0},
168  {1140,1.02643,0.34217,0.09158,0.00012,0,0,0},
169  {1145,1.02426,0.33869,0.08216,2e-05,0,0,0},
170  {1150,1.02085,0.33495,0.07213,-0.00011,0,0,0},
171  {1155,1.01581,0.33109,0.0622,-0.00026,0,0,0},
172  {1160,1.00873,0.32726,0.05308,-0.00044,0,0,0},
173  {1165,0.99942,0.32354,0.04523,-0.00065,0,0,0},
174  {1170,0.98823,0.31998,0.03857,-9e-04,0,0,0},
175  {1175,0.97557,0.31661,0.03292,-0.00117,0,0,0},
176  {1180,0.96183,0.31345,0.02813,-0.00148,0,0,0},
177  {1185,0.94743,0.31053,0.02403,-0.00183,0,0,0},
178  {1190,0.93279,0.30788,0.02045,-0.00221,1e-05,-1e-05,0},
179  {1195,0.91829,0.30552,0.01722,-0.00263,1e-05,-1e-05,0},
180  {1200,0.90437,0.30349,0.0142,-0.00309,1e-05,-1e-05,0},
181  {1205,0.89141,0.30181,0.0112,-0.0036,1e-05,-1e-05,0},
182  {1210,0.87983,0.30052,0.00806,-0.00414,1e-05,-1e-05,0},
183  {1215,0.87005,0.29962,0.00463,-0.00473,1e-05,-1e-05,0},
184  {1220,0.86245,0.29917,0.00072,-0.00537,1e-05,-1e-05,0},
185  {1225,0.85746,0.29917,-0.00381,-0.00605,1e-05,-1e-05,0},
186  {1230,0.85526,0.29965,-0.00909,-0.00678,1e-05,-1e-05,0},
187  {1235,0.85551,0.30056,-0.01516,-0.00754,1e-05,-2e-05,0},
188  {1240,0.85783,0.30188,-0.02202,-0.00833,1e-05,-2e-05,0},
189  {1245,0.86183,0.30357,-0.02969,-0.00912,0,-2e-05,0},
190  {1250,0.86711,0.30559,-0.03818,-0.00992,0,-2e-05,0},
191  {1255,0.87328,0.30789,-0.0475,-0.0107,0,-2e-05,0},
192  {1260,0.87995,0.31046,-0.05768,-0.01146,0,-2e-05,0},
193  {1265,0.88672,0.31324,-0.06872,-0.01218,0,-2e-05,0},
194  {1270,0.89321,0.3162,-0.08064,-0.01285,0,-2e-05,0},
195  {1275,0.89901,0.31931,-0.09346,-0.01347,0,-2e-05,0},
196  {1280,0.90375,0.32252,-0.10718,-0.01401,0,-2e-05,0},
197  {1285,0.90702,0.3258,-0.12182,-0.01447,-1e-05,-2e-05,0},
198  {1290,0.90844,0.32912,-0.13739,-0.01483,-1e-05,-2e-05,-1e-05},
199  {1295,0.90798,0.33251,-0.15381,-0.01508,-1e-05,-2e-05,-1e-05},
200  {1300,0.90606,0.33607,-0.17084,-0.01516,-2e-05,-2e-05,-1e-05},
201  {1305,0.90314,0.33995,-0.18827,-0.01504,-2e-05,-2e-05,-1e-05},
202  {1310,0.89965,0.34425,-0.20584,-0.01466,-2e-05,-2e-05,-1e-05},
203  {1315,0.89606,0.34911,-0.22335,-0.01398,-3e-05,-1e-05,-1e-05},
204  {1320,0.89281,0.35464,-0.24054,-0.01297,-3e-05,-1e-05,-2e-05},
205  {1325,0.89036,0.36097,-0.2572,-0.01157,-4e-05,-1e-05,-2e-05},
206  {1330,0.88915,0.36823,-0.27308,-0.00975,-5e-05,0,-2e-05},
207  {1335,0.88959,0.37652,-0.28798,-0.00745,-5e-05,0,-2e-05},
208  {1340,0.89162,0.38589,-0.30176,-0.00461,-6e-05,1e-05,-3e-05},
209  {1345,0.89499,0.39635,-0.31438,-0.00118,-6e-05,2e-05,-3e-05},
210  {1350,0.89948,0.40792,-0.32575,0.00291,-7e-05,3e-05,-3e-05},
211  {1355,0.90484,0.42061,-0.3358,0.00773,-8e-05,4e-05,-4e-05},
212  {1360,0.91082,0.43443,-0.34446,0.01335,-8e-05,5e-05,-4e-05},
213  {1365,0.91718,0.44939,-0.35167,0.01983,-9e-05,7e-05,-4e-05},
214  {1370,0.92368,0.46552,-0.35736,0.02723,-1e-04,8e-05,-4e-05},
215  {1375,0.93008,0.48281,-0.36144,0.03562,-1e-04,1e-04,-5e-05},
216  {1380,0.93614,0.50128,-0.36387,0.04506,-0.00011,0.00012,-5e-05},
217  {1385,0.94161,0.52096,-0.36455,0.05562,-0.00011,0.00014,-5e-05},
218  {1390,0.94624,0.54184,-0.36343,0.06737,-0.00012,0.00016,-5e-05},
219  {1395,0.94981,0.56394,-0.36043,0.08036,-0.00012,0.00018,-5e-05},
220  {1400,0.95212,0.58725,-0.35554,0.09464,-0.00012,0.00021,-5e-05},
221  {1405,0.95323,0.61167,-0.34898,0.11013,-0.00012,0.00023,-5e-05},
222  {1410,0.95325,0.63706,-0.34102,0.12669,-0.00012,0.00026,-5e-05},
223  {1415,0.95229,0.6633,-0.33195,0.14424,-0.00012,0.00029,-5e-05},
224  {1420,0.95046,0.69026,-0.32204,0.16264,-0.00011,0.00032,-4e-05},
225  {1425,0.94788,0.7178,-0.31156,0.1818,-1e-04,0.00035,-4e-05},
226  {1430,0.94466,0.74581,-0.30079,0.20158,-9e-05,0.00038,-3e-05},
227  {1435,0.9409,0.77415,-0.29002,0.22189,-8e-05,0.00041,-3e-05},
228  {1440,0.93673,0.80269,-0.2795,0.24261,-6e-05,0.00044,-2e-05},
229  {1445,0.93225,0.8313,-0.26953,0.26362,-4e-05,0.00047,-1e-05},
230  {1450,0.92757,0.85985,-0.26038,0.28481,-1e-05,5e-04,1e-05},
231  {1455,0.92282,0.88822,-0.25232,0.30607,2e-05,0.00052,2e-05},
232  {1460,0.91806,0.91629,-0.24556,0.3273,6e-05,0.00055,4e-05},
233  {1465,0.9133,0.94406,-0.24008,0.34847,1e-04,0.00057,6e-05},
234  {1470,0.90849,0.97153,-0.23576,0.36956,0.00015,0.00058,8e-05},
235  {1475,0.90359,0.99871,-0.2325,0.39056,0.00021,0.00059,0.00011},
236  {1480,0.89857,1.02561,-0.2302,0.41146,0.00027,0.00059,0.00014},
237  {1485,0.89338,1.05224,-0.22876,0.43223,0.00034,0.00059,0.00017},
238  {1490,0.88799,1.0786,-0.22808,0.45286,0.00042,0.00057,2e-04},
239  {1495,0.88235,1.1047,-0.22805,0.47333,5e-04,0.00054,0.00024},
240  {1500,0.87644,1.13056,-0.22858,0.49364,0.00059,0.00049,0.00028},
241  {1505,0.87021,1.15617,-0.22955,0.51376,0.00069,0.00043,0.00033},
242  {1510,0.86362,1.18155,-0.23086,0.53368,8e-04,0.00035,0.00037},
243  {1515,0.85663,1.20671,-0.23242,0.55338,0.00092,0.00026,0.00043},
244  {1520,0.84921,1.23163,-0.23416,0.57284,0.00104,0.00014,0.00048},
245  {1525,0.84136,1.25616,-0.23614,0.59195,0.00118,1e-05,0.00054},
246  {1530,0.83308,1.28016,-0.23847,0.6106,0.00132,-0.00015,6e-04},
247  {1535,0.82437,1.30346,-0.24125,0.62869,0.00146,-0.00034,0.00066},
248  {1540,0.81522,1.32589,-0.24459,0.6461,0.0016,-0.00055,0.00072},
249  {1545,0.80565,1.34729,-0.2486,0.66271,0.00175,-8e-04,0.00079},
250  {1550,0.79564,1.3675,-0.25336,0.67841,0.00189,-0.00107,0.00085},
251  {1555,0.78521,1.38636,-0.259,0.6931,0.00204,-0.00138,0.00091},
252  {1560,0.77436,1.4037,-0.26561,0.70666,0.00217,-0.00173,0.00096},
253  {1565,0.76308,1.41937,-0.2733,0.71898,0.0023,-0.00211,0.00101},
254  {1570,0.75138,1.4332,-0.28217,0.72994,0.00243,-0.00254,0.00106},
255  {1575,0.73925,1.44502,-0.29233,0.73943,0.00254,-0.00301,0.0011},
256  {1580,0.72671,1.45474,-0.30383,0.74739,0.00264,-0.00352,0.00114},
257  {1585,0.71375,1.46248,-0.31657,0.75388,0.00272,-0.00407,0.00116},
258  {1590,0.7004,1.46841,-0.33037,0.75901,0.00278,-0.00465,0.00118},
259  {1595,0.68665,1.47272,-0.3451,0.76289,0.0028,-0.00527,0.00118},
260  {1600,0.67251,1.4756,-0.36058,0.76565,0.00278,-0.00591,0.00115},
261  {1605,0.658,1.47722,-0.37667,0.76737,0.00272,-0.00657,0.00111},
262  {1610,0.64313,1.47777,-0.3932,0.76819,0.00259,-0.00725,0.00104},
263  {1615,0.6279,1.47743,-0.41002,0.76821,0.00241,-0.00795,0.00094},
264  {1620,0.61232,1.47637,-0.42696,0.76754,0.00215,-0.00865,0.00081},
265  {1625,0.59641,1.47479,-0.44387,0.76629,0.00182,-0.00935,0.00064},
266  {1630,0.58017,1.47287,-0.4606,0.76457,0.0014,-0.01005,0.00043},
267  {1635,0.56361,1.47078,-0.47698,0.7625,0.00089,-0.01075,0.00018},
268  {1640,0.54674,1.4687,-0.49286,0.76019,0.00028,-0.01143,-0.00011},
269  {1645,0.52957,1.46683,-0.50808,0.75774,-0.00044,-0.0121,-0.00045},
270  {1650,0.51216,1.46529,-0.5225,0.75525,-0.00127,-0.01274,-0.00085},
271  {1655,0.49484,1.46406,-0.53609,0.75274,-0.00222,-0.01333,-0.0013},
272  {1660,0.47796,1.46306,-0.54882,0.75023,-0.0033,-0.01384,-0.00181},
273  {1665,0.46188,1.4622,-0.56068,0.7477,-0.00452,-0.01424,-0.00237},
274  {1670,0.44698,1.46141,-0.57165,0.74517,-0.00587,-0.01449,-0.003},
275  {1675,0.43363,1.46061,-0.58171,0.74264,-0.00737,-0.01456,-0.00369},
276  {1680,0.42219,1.45972,-0.59084,0.74012,-0.00902,-0.01443,-0.00444},
277  {1685,0.41302,1.45866,-0.59903,0.73761,-0.01083,-0.01405,-0.00527},
278  {1690,0.40649,1.45737,-0.60625,0.73511,-0.01281,-0.01341,-0.00616},
279  {1695,0.40298,1.45575,-0.6125,0.73264,-0.01495,-0.01246,-0.00712},
280  {1700,0.40284,1.45373,-0.61774,0.73018,-0.01727,-0.01118,-0.00816},
281  {1705,0.40645,1.45123,-0.62196,0.72776,-0.01978,-0.00954,-0.00928},
282  {1710,0.41402,1.44822,-0.62516,0.72539,-0.02247,-0.00749,-0.01047},
283  {1715,0.42517,1.44489,-0.6274,0.72323,-0.02531,-0.005,-0.01172},
284  {1720,0.43937,1.44146,-0.62876,0.72143,-0.02827,-0.00201,-0.01302},
285  {1725,0.45609,1.43817,-0.62931,0.72016,-0.0313,0.00154,-0.01434},
286  {1730,0.4748,1.43525,-0.62913,0.71959,-0.03438,0.00569,-0.01566},
287  {1735,0.49497,1.43292,-0.6283,0.71989,-0.03746,0.01051,-0.01698},
288  {1740,0.51607,1.43142,-0.6269,0.72121,-0.04051,0.01603,-0.01826},
289  {1745,0.53758,1.43098,-0.62499,0.72373,-0.04349,0.02233,-0.0195},
290  {1750,0.55895,1.43183,-0.62267,0.7276,-0.04637,0.02944,-0.02067},
291  {1755,0.57967,1.4342,-0.61999,0.73301,-0.0491,0.03742,-0.02175},
292  {1760,0.5992,1.43833,-0.61705,0.74011,-0.05165,0.04632,-0.02273},
293  {1765,0.61701,1.44443,-0.61392,0.74907,-0.05398,0.0562,-0.0236},
294  {1770,0.63257,1.45275,-0.61067,0.76005,-0.05606,0.06712,-0.02432},
295  {1775,0.64535,1.46351,-0.60738,0.77322,-0.05785,0.07911,-0.02489},
296  {1780,0.65496,1.47687,-0.60408,0.78868,-0.05929,0.09221,-0.02527},
297  {1785,0.66162,1.49263,-0.60057,0.80621,-0.06023,0.10632,-0.02541},
298  {1790,0.66566,1.51054,-0.5966,0.82552,-0.06051,0.12132,-0.02522},
299  {1795,0.66742,1.53032,-0.59194,0.84633,-0.05997,0.13708,-0.02464},
300  {1800,0.66727,1.55172,-0.58634,0.86836,-0.05843,0.15349,-0.02359},
301  {1805,0.66553,1.57446,-0.57955,0.89131,-0.05573,0.17042,-0.02199},
302  {1810,0.66256,1.59828,-0.57132,0.9149,-0.05171,0.18775,-0.01978},
303  {1815,0.6587,1.62291,-0.56141,0.93885,-0.04619,0.20535,-0.01688},
304  {1820,0.6543,1.64809,-0.54958,0.96288,-0.03901,0.2231,-0.01321},
305  {1825,0.64971,1.67356,-0.53558,0.98668,-0.03001,0.24088,-0.00871},
306  {1830,0.64526,1.69904,-0.51916,1.00999,-0.01902,0.25856,-0.00329},
307  {1835,0.64131,1.72427,-0.50008,1.0325,-0.00587,0.27603,0.00312},
308  {1840,0.63821,1.74898,-0.47809,1.05395,0.00961,0.29315,0.01058},
309  {1845,0.63629,1.77291,-0.45295,1.07404,0.02757,0.30981,0.01917},
310  {1850,0.63591,1.79579,-0.42441,1.09248,0.04819,0.32589,0.02898},
311  {1855,0.63741,1.81736,-0.39222,1.10899,0.07164,0.34125,0.04006},
312  {1860,0.64113,1.83735,-0.35615,1.12328,0.09807,0.35578,0.0525},
313  {1865,0.64742,1.8555,-0.31594,1.13507,0.12766,0.36936,0.06637},
314  {1870,0.65664,1.87153,-0.27136,1.14408,0.16058,0.38186,0.08174},
315  {1875,0.66911,1.88519,-0.22214,1.15001,0.19698,0.39315,0.09868},
316  {1880,0.68509,1.89627,-0.16819,1.15266,0.23695,0.40315,0.11724},
317  {1885,0.70432,1.90488,-0.10992,1.15218,0.28023,0.4119,0.1373},
318  {1890,0.72647,1.91118,-0.04788,1.14879,0.3265,0.41946,0.15872},
319  {1895,0.75118,1.91533,0.01738,1.14272,0.37541,0.42591,0.18135},
320  {1900,0.77812,1.91751,0.08531,1.13419,0.42662,0.43131,0.20505},
321  {1905,0.80692,1.91788,0.15535,1.12342,0.47979,0.43573,0.22967},
322  {1910,0.83725,1.91661,0.22697,1.11064,0.53459,0.43925,0.25507},
323  {1915,0.86875,1.91386,0.2996,1.09608,0.59068,0.44194,0.28109},
324  {1920,0.90109,1.90981,0.37269,1.07996,0.64773,0.44385,0.30761},
325  {1925,0.9339,1.90461,0.4457,1.06251,0.70538,0.44508,0.33446},
326  {1930,0.96685,1.89844,0.51808,1.04394,0.76331,0.44567,0.36151},
327  {1935,0.99958,1.89146,0.58927,1.02449,0.82118,0.44571,0.38861},
328  {1940,1.03175,1.88384,0.65873,1.00439,0.87864,0.44526,0.41561},
329  {1945,1.06301,1.87575,0.72589,0.98385,0.93537,0.44439,0.44237},
330  {1950,1.09301,1.86735,0.79022,0.9631,0.99102,0.44318,0.46874},
331  {1955,1.1214,1.85881,0.85116,0.94236,1.04526,0.44169,0.49458},
332  {1960,1.14785,1.8503,0.90817,0.92187,1.09774,0.43999,0.51975},
333  {1965,1.17199,1.84198,0.96068,0.90184,1.14813,0.43816,0.54408},
334  {1970,1.19348,1.83402,1.00815,0.88251,1.1961,0.43625,0.56745},
335  {1975,1.21198,1.82659,1.05004,0.86409,1.2413,0.43435,0.58971},
336  {1980,1.22713,1.81986,1.08578,0.84681,1.28339,0.43252,0.6107},
337  {1985,1.23859,1.81398,1.11483,0.8309,1.32204,0.43084,0.63029},
338  {1990,1.24601,1.80914,1.13664,0.81658,1.35691,0.42936,0.64833},
339  {1995,1.24905,1.80549,1.15065,0.80408,1.38766,0.42817,0.66467},
340  {2000,1.24735,1.8032,1.15632,0.79362,1.41396,0.42732,0.67917},
341  {2005,1.2407,1.80242,1.15334,0.78541,1.43554,0.42692,0.6917},
342  {2010,1.22949,1.80321,1.14236,0.77959,1.45242,0.42713,0.70224},
343  {2015,1.21419,1.80564,1.12429,0.7763,1.4647,0.42812,0.71076},
344  {2020,1.19533,1.80973,1.10003,0.77565,1.47248,0.4301,0.71727},
345  {2025,1.17341,1.81555,1.07046,0.77778,1.47586,0.43325,0.72176},
346  {2030,1.14892,1.82314,1.0365,0.78281,1.47493,0.43775,0.7242},
347  {2035,1.12238,1.83254,0.99904,0.79087,1.4698,0.44379,0.7246},
348  {2040,1.09429,1.84381,0.95898,0.8021,1.46056,0.45155,0.72294},
349  {2045,1.06515,1.85699,0.91723,0.8166,1.44731,0.46123,0.71921},
350  {2050,1.03547,1.87214,0.87467,0.83452,1.43015,0.473,0.71341},
351  {2055,1.00576,1.88929,0.83222,0.85598,1.40917,0.48707,0.70551},
352  {2060,0.97644,1.90846,0.79062,0.88103,1.38453,0.50353,0.69555},
353  {2065,0.94767,1.92947,0.75007,0.90941,1.35653,0.52226,0.6837},
354  {2070,0.91954,1.95211,0.71063,0.94081,1.32555,0.54305,0.67014},
355  {2075,0.89213,1.97615,0.67234,0.97489,1.29194,0.5657,0.65508},
356  {2080,0.86552,2.0014,0.63526,1.01133,1.25607,0.58999,0.63872},
357  {2085,0.83981,2.02762,0.59943,1.04978,1.21831,0.61572,0.62125},
358  {2090,0.81506,2.0546,0.56491,1.08994,1.17903,0.64269,0.60288},
359  {2095,0.79138,2.08213,0.53175,1.13146,1.13857,0.67069,0.5838},
360  {2100,0.76883,2.11,0.5,1.17402,1.09732,0.69951,0.56422},
361  {2105,0.74751,2.13798,0.46972,1.2173,1.05563,0.72895,0.54432},
362  {2110,0.72751,2.16586,0.44094,1.26096,1.01387,0.7588,0.52432},
363  {2115,0.70889,2.19343,0.41374,1.30467,0.9724,0.78886,0.5044},
364  {2120,0.69174,2.22049,0.38814,1.34815,0.93155,0.81894,0.48475},
365  {2125,0.67605,2.24701,0.36417,1.39128,0.89144,0.84898,0.46544},
366  {2130,0.6618,2.27294,0.34183,1.43399,0.85218,0.87891,0.44653},
367  {2135,0.64897,2.29827,0.32114,1.47621,0.81386,0.90871,0.42808},
368  {2140,0.63754,2.32297,0.3021,1.51786,0.77658,0.9383,0.41014},
369  {2145,0.62749,2.34703,0.28473,1.55887,0.74042,0.96766,0.39278},
370  {2150,0.61881,2.37041,0.26902,1.59918,0.70547,0.99672,0.37603},
371  {2155,0.61148,2.39309,0.25499,1.63869,0.67185,1.02545,0.35997},
372  {2160,0.60547,2.41504,0.24265,1.67735,0.63963,1.05379,0.34465},
373  {2165,0.60077,2.43625,0.23201,1.71507,0.60891,1.08169,0.33013},
374  {2170,0.59735,2.45669,0.22308,1.75179,0.57979,1.10911,0.31645},
375  {2175,0.59521,2.47634,0.21586,1.78743,0.55236,1.136,0.30368},
376  {2180,0.59431,2.49517,0.21037,1.82192,0.52671,1.1623,0.29188},
377  {2185,0.59465,2.51315,0.20661,1.85518,0.50293,1.18798,0.2811},
378  {2190,0.59619,2.53027,0.20459,1.88715,0.48113,1.21298,0.27139},
379  {2195,0.59893,2.5465,0.20433,1.91774,0.46139,1.23725,0.26282},
380  {2200,0.60284,2.56181,0.20582,1.9469,0.44381,1.26075,0.25544},
381  {2205,0.60791,2.57618,0.20909,1.97453,0.42849,1.28343,0.2493},
382  {2210,0.61411,2.5896,0.21413,2.00058,0.41551,1.30524,0.24447},
383  {2215,0.62143,2.60202,0.22096,2.02496,0.40496,1.32612,0.24099},
384  {2220,0.62984,2.61344,0.22959,2.04761,0.39696,1.34604,0.23893},
385  {2225,0.63934,2.62382,0.24002,2.06845,0.39157,1.36495,0.23834},
386  {2230,0.64989,2.63314,0.25227,2.0874,0.38891,1.38278,0.23928},
387  {2235,0.66148,2.64138,0.26634,2.1044,0.38907,1.39951,0.24181},
388  {2240,0.67408,2.64852,0.28223,2.1194,0.3921,1.41509,0.24595},
389  {2245,0.68768,2.6546,0.29989,2.13244,0.39794,1.42955,0.25168},
390  {2250,0.70222,2.65965,0.31925,2.14359,0.40647,1.44294,0.25893},
391  {2255,0.71767,2.66371,0.34024,2.15294,0.4176,1.45528,0.26765},
392  {2260,0.73399,2.66683,0.36281,2.16054,0.43123,1.46664,0.27777},
393  {2265,0.75115,2.66903,0.38688,2.16649,0.44724,1.47703,0.28924},
394  {2270,0.76912,2.67037,0.41239,2.17083,0.46554,1.48652,0.30199},
395  {2275,0.78785,2.67088,0.43927,2.17366,0.48602,1.49514,0.31598},
396  {2280,0.80731,2.67059,0.46746,2.17504,0.50857,1.50293,0.33114},
397  {2285,0.82746,2.66956,0.49689,2.17504,0.5331,1.50993,0.34741},
398  {2290,0.84826,2.66782,0.5275,2.17375,0.55949,1.51618,0.36473},
399  {2295,0.86969,2.66541,0.55922,2.17122,0.58765,1.52173,0.38305},
400  {2300,0.89169,2.66236,0.59199,2.16753,0.61747,1.52662,0.40231},
401  {2305,0.91425,2.65872,0.62574,2.16276,0.64884,1.53089,0.42244},
402  {2310,0.93731,2.65453,0.6604,2.15697,0.68167,1.53458,0.44339},
403  {2315,0.96085,2.64983,0.69591,2.15025,0.71584,1.53773,0.46511},
404  {2320,0.98482,2.64466,0.7322,2.14266,0.75126,1.54038,0.48752},
405  {2325,1.00919,2.63905,0.76922,2.13427,0.78782,1.54258,0.51058},
406  {2330,1.03393,2.63304,0.80688,2.12516,0.82541,1.54437,0.53422},
407  {2335,1.05899,2.62669,0.84513,2.11541,0.86393,1.54578,0.55838},
408  {2340,1.08434,2.62001,0.88391,2.10507,0.90328,1.54686,0.58302},
409  {2345,1.10995,2.61307,0.92313,2.09423,0.94336,1.54766,0.60806},
410  {2350,1.13577,2.60588,0.96275,2.08296,0.98405,1.54821,0.63345},
411  {2355,1.16178,2.5985,1.0027,2.07133,1.02526,1.54855,0.65913},
412  {2360,1.18793,2.59097,1.0429,2.05941,1.06688,1.54873,0.68504},
413  {2365,1.21419,2.58332,1.0833,2.04728,1.10881,1.54878,0.71113},
414  {2370,1.24052,2.57559,1.12382,2.035,1.15094,1.54876,0.73733},
415  {2375,1.26688,2.56782,1.16441,2.02266,1.19317,1.54869,0.76358}};
416 
417  const G4int coef_ener = G4int((pLab-Legendre_coef[0][0])/5);
418  const G4double sup_ener = pLab/5. - coef_ener -Legendre_coef[0][0]/5;
419  //std::cout << "sup_ener\t" << sup_ener << std::endl;
420 
421 // assert(pLab >= Legendre_coef[coef_ener][0] && pLab < Legendre_coef[coef_ener+1][0]);
422 
423  // Legendre coefficient normalized
424  const G4double A0 = 1.;
425  const G4double A1 = (1-sup_ener)*Legendre_coef[coef_ener][1] + sup_ener*Legendre_coef[coef_ener+1][1];
426  const G4double A2 = (1-sup_ener)*Legendre_coef[coef_ener][2] + sup_ener*Legendre_coef[coef_ener+1][2];
427  const G4double A3 = (1-sup_ener)*Legendre_coef[coef_ener][3] + sup_ener*Legendre_coef[coef_ener+1][3];
428  const G4double A4 = (1-sup_ener)*Legendre_coef[coef_ener][4] + sup_ener*Legendre_coef[coef_ener+1][4];
429  const G4double A5 = (1-sup_ener)*Legendre_coef[coef_ener][5] + sup_ener*Legendre_coef[coef_ener+1][5];
430  const G4double A6 = (1-sup_ener)*Legendre_coef[coef_ener][6] + sup_ener*Legendre_coef[coef_ener+1][6];
431  const G4double A7 = (1-sup_ener)*Legendre_coef[coef_ener][7] + sup_ener*Legendre_coef[coef_ener+1][7];
432 
433  // Theoritical max if all Ai > 0 (often the case)
434  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);
435 
436  G4bool success = false;
437  G4int maxloop = 0;
438 
439  while(!success && maxloop < 1000){
440 
441  cos_theta = Random::shoot()*2-1.; // not optimized
442 
443  // Legendre Polynomial
444  G4double P0 = A0;
445  G4double P1 = A1*cos_theta;
446  G4double P2 = A2/2.*(3*std::pow(cos_theta,2)-1);
447  G4double P3 = A3/2.*(5*std::pow(cos_theta,3)-3*cos_theta);
448  G4double P4 = A4/8.*(35*std::pow(cos_theta,4)-30*std::pow(cos_theta,2)+3);
449  G4double P5 = A5/8.*(63*std::pow(cos_theta,5)-70*std::pow(cos_theta,3)+15*cos_theta);
450  G4double P6 = A6/16.*(231*std::pow(cos_theta,6)-315*std::pow(cos_theta,4)+105*std::pow(cos_theta,2)-5);
451  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);
452 
453  G4double P = (P0 + P1 + P2 + P3 + P4 + P5 + P6 + P7);
454 
455  if(Random::shoot()*A < P) success = true;
456  maxloop +=1 ;
457  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
458  }
459  sin_theta = std::sqrt(1-cos_theta*cos_theta);
460  }
461 
462  if(rho == 0) return ThreeVector(sin_theta*cos_phi,sin_theta*sin_phi,cos_theta);
463  // Rotation in the direction of the incident pi
464  const G4double px = x/r*cos_theta - y/rho*sin_theta*cos_phi + z/r*x/rho*sin_theta*sin_phi;
465  const G4double py = y/r*cos_theta + x/rho*sin_theta*cos_phi + z/r*y/rho*sin_theta*sin_phi;
466  const G4double pz = z/r*cos_theta - rho/r*sin_theta*sin_phi;
467 
468 
469  return ThreeVector(px,py,pz);
470  }
471 
472 }