ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLPhaseSpaceRauboldLynch.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLPhaseSpaceRauboldLynch.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 
39 #include "G4INCLRandom.hh"
40 #include "G4INCLKinematicsUtils.hh"
41 #include <algorithm>
42 #include <numeric>
43 #include <functional>
44 
45 namespace G4INCL {
46 
48  0.01,
49  0.017433288222,
50  0.0303919538231,
51  0.0529831690628,
52  0.0923670857187,
53  0.161026202756,
54  0.280721620394,
55  0.489390091848,
56  0.853167852417,
57  1.48735210729,
58  2.5929437974,
59  4.52035365636,
60  7.88046281567,
61  13.7382379588,
62  23.9502661999,
63  41.7531893656,
64  72.7895384398,
65  126.896100317,
66  221.221629107,
67  385.662042116,
68  672.33575365,
69  1172.10229753,
70  2043.35971786,
71  3562.24789026,
72  6210.16941892,
73  10826.3673387,
74  18873.9182214,
75  32903.4456231,
76  57361.5251045,
77  100000.0
78  };
79 
81  -4.69180212056,
82  -4.13600582212,
83  -3.58020940928,
84  -3.02441303018,
85  -2.46861663471,
86  -1.91282025644,
87  -1.35702379603,
88  -0.801227342882,
89  -0.245430866403,
90  0.310365269464,
91  0.866161720461,
92  1.42195839972,
93  1.97775459763,
94  2.5335510251,
95  3.08934734235,
96  3.64514378357,
97  4.20094013304,
98  4.75673663205,
99  5.3125329676,
100  5.86832950014,
101  6.42412601468,
102  6.97992197797,
103  7.53571856765,
104  8.09151503031,
105  8.64731144398,
106  9.20310774148,
107  9.7589041009,
108  10.314700625,
109  10.8704971207,
110  11.4262935304
111  };
112 
114  8.77396538453e-07,
115  1.52959067398e-06,
116  2.66657950812e-06,
117  4.6487249132e-06,
118  8.10425612766e-06,
119  1.41283832898e-05,
120  2.46304178003e-05,
121  4.2938917254e-05,
122  7.4856652043e-05,
123  0.00013049975904,
124  0.000227503991225,
125  0.000396614265067,
126  0.000691429079588,
127  0.00120538824295,
128  0.00210138806588,
129  0.00366341038188,
130  0.00638652890627,
131  0.0111338199161,
132  0.0194099091609,
133  0.0338378540766,
134  0.0589905062931,
135  0.102839849857,
136  0.179283674326,
137  0.312550396803,
138  0.544878115136,
139  0.949901722703,
140  1.65599105145,
141  2.88693692929,
142  5.03288035671,
143  8.77396538453
144  };
146  7.28682764024,
147  7.0089298602,
148  6.7310326046,
149  6.45313436085,
150  6.17523713298,
151  5.89734080335,
152  5.61944580818,
153  5.34155326611,
154  5.06366530628,
155  4.78578514804,
156  4.50791742491,
157  4.23007259554,
158  3.97566018117,
159  3.72116551935,
160  3.44355099732,
161  3.1746329047,
162  2.89761210229,
163  2.62143335535,
164  2.34649707964,
165  2.07366126445,
166  1.8043078447,
167  1.54056992953,
168  1.27913996411,
169  0.981358535322,
170  0.73694629682,
171  0.587421056631,
172  0.417178268911,
173  0.280881750176,
174  0.187480311508,
175  0.116321254846
176  };
177 
179 
181  nParticles(0),
182  sqrtS(0.),
183  availableEnergy(0.),
184  maxGeneratedWeight(0.)
185  {
186  std::vector<G4double> wMaxMasslessXV(wMaxMasslessX, wMaxMasslessX + wMaxNE);
187  std::vector<G4double> wMaxMasslessYV(wMaxMasslessY, wMaxMasslessY + wMaxNE);
188  wMaxMassless = new InterpolationTable(wMaxMasslessXV, wMaxMasslessYV);
189  std::vector<G4double> wMaxCorrectionXV(wMaxCorrectionX, wMaxCorrectionX + wMaxNE);
190  std::vector<G4double> wMaxCorrectionYV(wMaxCorrectionY, wMaxCorrectionY + wMaxNE);
191  wMaxCorrection = new InterpolationTable(wMaxCorrectionXV, wMaxCorrectionYV);
192 
193  // initialize the precalculated coefficients
194  prelog[0] = 0.;
195  for(size_t i=1; i<wMaxNP; ++i) {
196  prelog[i] = -std::log(G4double(i));
197  }
198  }
199 
201  delete wMaxMassless;
202  delete wMaxCorrection;
203  }
204 
206  maxGeneratedWeight = 0.;
207 
208  sqrtS = sqs;
209 
210  // initialize the structures containing the particle masses
211  initialize(particles);
212 
213  // initialize the maximum weight
214  const G4double weightMax = computeMaximumWeightParam();
215 
216  const G4int maxIter = 500;
217  G4int iter = 0;
218  G4double weight, r;
219  do {
220  weight = computeWeight();
222 // assert(weight<=weightMax);
223  r = Random::shoot();
224  } while(++iter<maxIter && r*weightMax>weight); /* Loop checking, 10.07.2015, D.Mancusi */
225 
226 #ifndef INCLXX_IN_GEANT4_MODE
227  if(iter>=maxIter) {
228  INCL_WARN("Number of tries exceeded in PhaseSpaceRauboldLynch::generate()\n"
229  << "nParticles=" << nParticles << ", sqrtS=" << sqrtS << '\n');
230  }
231 #endif
232 
233  generateEvent(particles);
234  }
235 
237  nParticles = particles.size();
238 // assert(nParticles>2);
239 
240  // masses and sum of masses
241  masses.resize(nParticles);
242  sumMasses.resize(nParticles);
243  std::transform(particles.begin(), particles.end(), masses.begin(), std::mem_fn(&Particle::getMass));
244  std::partial_sum(masses.begin(), masses.end(), sumMasses.begin());
245 
246  // sanity check
248 // assert(availableEnergy>-1.e-5);
249  if(availableEnergy<0.)
250  availableEnergy = 0.;
251 
252  rnd.resize(nParticles);
253  invariantMasses.resize(nParticles);
254  momentaCM.resize(nParticles-1);
255  }
256 
258  // compute the maximum weight
259  G4double eMMax = sqrtS + masses[0];
260  G4double eMMin = 0.;
261  G4double wMax = 1.;
262  for(size_t i=1; i<nParticles; i++) {
263  eMMin += masses[i-1];
264  eMMax += masses[i];
265  wMax *= KinematicsUtils::momentumInCM(eMMax, eMMin, masses[i]);
266  }
267  return wMax;
268  }
269 
271 #ifndef INCLXX_IN_GEANT4_MODE
272  if(nParticles>=wMaxNP) {
273  INCL_WARN("The requested number of particles (" << nParticles << ") requires extrapolation the tables in PhaseSpaceRauboldLynch. YMMV." << '\n');
274  }
275 
277  INCL_WARN("The requested available energy (" << availableEnergy << " MeV) requires extrapolation the tables in PhaseSpaceRauboldLynch. YMMV." << '\n');
278  }
279 #endif
280 
281  // compute the maximum weight
282  const G4double logMassless = ((*wMaxMassless)(availableEnergy)+prelog[nParticles])*(nParticles-1);
283  const G4double reducedSqrtS = availableEnergy/sumMasses[nParticles-1];
284  const G4double correction = (*wMaxCorrection)(reducedSqrtS);
285  const G4double wMax = std::exp(correction*G4double(nParticles-1) + logMassless + wMaxInterpolationMargin);
286  if(wMax>0.)
287  return wMax;
288  else
289  return computeMaximumWeightNaive();
290  }
291 
293  // Generate nParticles-2 sorted random numbers, bracketed by 0 and 1
294  rnd[0] = 0.;
295  for(size_t i=1; i<nParticles-1; ++i)
296  rnd[i] = Random::shoot();
297  rnd[nParticles-1] = 1.;
298  std::sort(rnd.begin()+1, rnd.begin()+nParticles-1);
299 
300  // invariant masses
301  for(size_t i=0; i<nParticles; ++i)
303 
304  // compute the CM momenta and the weight for the current event
306  momentaCM[0] = weight;
307  for(size_t i=1; i<nParticles-1; ++i) {
308  G4double momentumCM;
309 // assert(invariantMasses[i+1]-invariantMasses[i]-masses[i+1]> - 1.E-8);
310  if(invariantMasses[i+1]-invariantMasses[i]-masses[i+1] < 0.) momentumCM = 0.;
311  else momentumCM = KinematicsUtils::momentumInCM(invariantMasses[i+1], invariantMasses[i], masses[i+1]);
312  momentaCM[i] = momentumCM;
313  weight *= momentumCM;
314  }
315 
316  return weight;
317  }
318 
320  Particle *p = particles[0];
322  p->setMomentum(momentum);
324 
325  ThreeVector boostV;
326 
327  for(size_t i=1; i<nParticles; ++i) {
328  p = particles[i];
329  p->setMomentum(-momentum);
331 
332  if(i==nParticles-1)
333  break;
334 
335  momentum = Random::normVector(momentaCM[i]);
336 
337  const G4double iM = invariantMasses[i];
338  const G4double recoilE = std::sqrt(momentum.mag2() + iM*iM);
339  boostV = -momentum/recoilE;
340  for(size_t j=0; j<=i; ++j)
341  particles[j]->boost(boostV);
342  }
343  }
344 
346  return maxGeneratedWeight;
347  }
348 }