Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4CascadeFinalStateAlgorithm.cc
Go to the documentation of this file.
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 // $Id$
27 // Author: Michael Kelsey (SLAC)
28 // Date: 15 April 2013
29 //
30 // Description: Subclass of models/util G4VHadDecayAlgorithm which uses
31 // old INUCL parametrizations for momentum and angular
32 // distributions.
33 //
34 // 20130509 BUG FIX: Two-body momentum vector should be rotated into
35 // collision axis; three-body "final" vector needs to be rotated
36 // into axis of rest of system. Tweak some diagnostic messages
37 // to match old G4EPCollider version.
38 // 20130612 BUG FIX: Create separate FillDirThreeBody(), FillDirManyBody()
39 // in order to reporoduce old method: N-body states are filled
40 // from first to last, while three-body starts with the last.
41 // 20130702 M. Kelsey: Copy phase-space algorithm from Kopylov; use if
42 // runtime envvar G4CASCADE_USE_PHASESPACE is set
43 // 20140627 BUG FIX: Use ".c_str()" in diagnostics to avoid IBM XL error.
44 // 20150608 M. Kelsey -- Label all while loops as terminating.
45 // 20150619 M. Kelsey -- Replace std::exp with G4Exp
46 
48 #include "G4CascadeParameters.hh"
49 #include "G4Exp.hh"
52 #include "G4LorentzConvertor.hh"
54 #include "G4Pow.hh"
55 #include "G4TwoBodyAngularDist.hh"
56 #include "G4VMultiBodyMomDst.hh"
57 #include "G4VTwoBodyAngDst.hh"
58 #include "Randomize.hh"
59 #include <vector>
60 #include <numeric>
61 #include <cmath>
62 
63 using namespace G4InuclSpecialFunctions;
64 
65 
66 // Cut-offs and iteration limits for generation
67 
68 const G4double G4CascadeFinalStateAlgorithm::maxCosTheta = 0.9999;
69 const G4double G4CascadeFinalStateAlgorithm::oneOverE = 0.3678794;
70 const G4double G4CascadeFinalStateAlgorithm::small = 1.e-10;
71 const G4int G4CascadeFinalStateAlgorithm::itry_max = 10;
72 
73 
74 // Constructor and destructor
75 
77  : G4VHadDecayAlgorithm("G4CascadeFinalStateAlgorithm"),
78  momDist(0), angDist(0), multiplicity(0), bullet_ekin(0.) {;}
79 
81 
86  toSCM.setVerbose(verbose);
87 }
88 
89 
90 // Select distributions to be used for next interaction
91 
95  const std::vector<G4int>& particle_kinds) {
96  if (GetVerboseLevel()>1)
97  G4cout << " >>> " << GetName() << "::Configure" << G4endl;
98 
99  // Identify initial and final state (if two-body) for algorithm selection
100  multiplicity = particle_kinds.size();
101  G4int is = bullet->type() * target->type();
102  G4int fs = (multiplicity==2) ? particle_kinds[0]*particle_kinds[1] : 0;
103 
104  ChooseGenerators(is, fs);
105 
106  // Save kinematics for use with distributions
107  SaveKinematics(bullet, target);
108 
109  // Save particle types for use with distributions
110  kinds = particle_kinds;
111 }
112 
113 // Save kinematics for use with generators
114 
118  if (GetVerboseLevel()>1)
119  G4cout << " >>> " << GetName() << "::SaveKinematics" << G4endl;
120 
121  if (target->nucleon()) { // Which particle originated in nucleus?
122  toSCM.setBullet(bullet);
123  toSCM.setTarget(target);
124  } else {
125  toSCM.setBullet(target);
126  toSCM.setTarget(bullet);
127  }
128 
129  toSCM.toTheCenterOfMass();
130 
131  bullet_ekin = toSCM.getKinEnergyInTheTRS();
132 }
133 
134 
135 // Select generator based on initial and final state
136 
138  if (GetVerboseLevel()>1)
139  G4cout << " >>> " << GetName() << "::ChooseGenerators"
140  << " is " << is << " fs " << fs << G4endl;
141 
142  // Get generators for momentum and angle
143  if (G4CascadeParameters::usePhaseSpace()) momDist = 0;
144  else momDist = G4MultiBodyMomentumDist::GetDist(is, multiplicity);
145 
146  if (fs > 0 && multiplicity == 2) {
147  G4int kw = (fs==is) ? 1 : 2;
148  angDist = G4TwoBodyAngularDist::GetDist(is, fs, kw);
149  } else if (multiplicity == 3) {
150  angDist = G4TwoBodyAngularDist::GetDist(is);
151  } else {
152  angDist = 0;
153  }
154 
155  if (GetVerboseLevel()>1) {
156  G4cout << " " << (momDist?momDist->GetName().c_str():"") << " "
157  << (angDist?angDist->GetName().c_str():"") << G4endl;
158  }
159 }
160 
161 
162 // Two-body generation uses angular-distribution function
163 
165 GenerateTwoBody(G4double initialMass, const std::vector<G4double>& masses,
166  std::vector<G4LorentzVector>& finalState) {
167  if (GetVerboseLevel()>1)
168  G4cout << " >>> " << GetName() << "::GenerateTwoBody" << G4endl;
169 
170  finalState.clear(); // Initialization and sanity checks
171 
172  if (multiplicity != 2) return;
173 
174  // Generate momentum vector in CMS for back-to-back particles
175  G4double pscm = TwoBodyMomentum(initialMass, masses[0], masses[1]);
176 
177  G4double costh = angDist ? angDist->GetCosTheta(bullet_ekin, pscm)
178  : (2.*G4UniformRand() - 1.);
179 
180  mom.setRThetaPhi(pscm, std::acos(costh), UniformPhi());
181 
182  if (GetVerboseLevel()>3) { // Copied from old G4EPCollider
183  G4cout << " Particle kinds = " << kinds[0] << " , " << kinds[1]
184  << "\n pmod " << pscm
185  << "\n before rotation px " << mom.x() << " py " << mom.y()
186  << " pz " << mom.z() << G4endl;
187  }
188 
189  finalState.resize(2); // Allows filling by index
190 
191  finalState[0].setVectM(mom, masses[0]);
192  finalState[0] = toSCM.rotate(finalState[0]);
193 
194  if (GetVerboseLevel()>3) { // Copied from old G4EPCollider
195  G4cout << " after rotation px " << finalState[0].x() << " py "
196  << finalState[0].y() << " pz " << finalState[0].z() << G4endl;
197  }
198 
199  finalState[1].setVectM(-finalState[0].vect(), masses[1]);
200 }
201 
202 
203 // N-body generation uses momentum-modulus distribution, computed angles
204 
206 GenerateMultiBody(G4double initialMass, const std::vector<G4double>& masses,
207  std::vector<G4LorentzVector>& finalState) {
208  if (GetVerboseLevel()>1)
209  G4cout << " >>> " << GetName() << "::GenerateMultiBody" << G4endl;
210 
212  FillUsingKopylov(initialMass, masses, finalState);
213  return;
214  }
215 
216  finalState.clear(); // Initialization and sanity checks
217  if (multiplicity < 3) return;
218  if (!momDist) return;
219 
220  G4int itry = -1; /* Loop checking 08.06.2015 MHK */
221  while ((G4int)finalState.size() != multiplicity && ++itry < itry_max) {
222  FillMagnitudes(initialMass, masses);
223  FillDirections(initialMass, masses, finalState);
224  }
225 }
226 
227 
229 FillMagnitudes(G4double initialMass, const std::vector<G4double>& masses) {
230  if (GetVerboseLevel()>1)
231  G4cout << " >>> " << GetName() << "::FillMagnitudes" << G4endl;
232 
233  modules.clear(); // Initialization and sanity checks
234  if (!momDist) return;
235 
236  modules.resize(multiplicity,0.); // Pre-allocate to avoid resizing
237 
238  G4double mass_last = masses.back();
239  G4double pmod = 0.;
240 
241  if (GetVerboseLevel() > 3){
242  G4cout << " knd_last " << kinds.back() << " mass_last "
243  << mass_last << G4endl;
244  }
245 
246  G4int itry = -1;
247  while (++itry < itry_max) { /* Loop checking 08.06.2015 MHK */
248  if (GetVerboseLevel() > 3) {
249  G4cout << " itry in fillMagnitudes " << itry << G4endl;
250  }
251 
252  G4double eleft = initialMass;
253 
254  G4int i; // For access outside of loop
255  for (i=0; i < multiplicity-1; i++) {
256  pmod = momDist->GetMomentum(kinds[i], bullet_ekin);
257 
258  if (pmod < small) break;
259  eleft -= std::sqrt(pmod*pmod + masses[i]*masses[i]);
260 
261  if (GetVerboseLevel() > 3) {
262  G4cout << " kp " << kinds[i] << " pmod " << pmod
263  << " mass2 " << masses[i]*masses[i] << " eleft " << eleft
264  << "\n x1 " << eleft - mass_last << G4endl;
265  }
266 
267  if (eleft <= mass_last) break;
268 
269  modules[i] = pmod;
270  }
271 
272  if (i < multiplicity-1) continue; // Failed to generate full kinematics
273 
274  G4double plast = eleft * eleft - masses.back()*masses.back();
275  if (GetVerboseLevel() > 2) G4cout << " plast ** 2 " << plast << G4endl;
276 
277  if (plast <= small) continue; // Not enough momentum left over
278 
279  plast = std::sqrt(plast); // Final momentum is what's left over
280  modules.back() = plast;
281 
282  if (multiplicity > 3 || satisfyTriangle(modules)) break; // Successful
283  } // while (itry < itry_max)
284 
285  if (itry >= itry_max) { // Too many attempts
286  if (GetVerboseLevel() > 2)
287  G4cerr << " Unable to generate momenta for multiplicity "
288  << multiplicity << G4endl;
289 
290  modules.clear(); // Something went wrong, throw away partial
291  }
292 }
293 
294 // For three-body states, check kinematics of momentum moduli
295 
297 satisfyTriangle(const std::vector<G4double>& pmod) const {
298  if (GetVerboseLevel() > 3)
299  G4cout << " >>> " << GetName() << "::satisfyTriangle" << G4endl;
300 
301  return ( (pmod.size() != 3) ||
302  !(pmod[0] < std::fabs(pmod[1] - pmod[2]) ||
303  pmod[0] > pmod[1] + pmod[2] ||
304  pmod[1] < std::fabs(pmod[0] - pmod[2]) ||
305  pmod[1] > pmod[0] + pmod[2] ||
306  pmod[2] < std::fabs(pmod[0] - pmod[1]) ||
307  pmod[2] > pmod[1] + pmod[0])
308  );
309 }
310 
311 // Generate momentum directions into final state
312 
314 FillDirections(G4double initialMass, const std::vector<G4double>& masses,
315  std::vector<G4LorentzVector>& finalState) {
316  if (GetVerboseLevel()>1)
317  G4cout << " >>> " << GetName() << "::FillDirections" << G4endl;
318 
319  finalState.clear(); // Initialization and sanity check
320  if ((G4int)modules.size() != multiplicity) return;
321 
322  // Different order of processing for three vs. N body
323  if (multiplicity == 3)
324  FillDirThreeBody(initialMass, masses, finalState);
325  else
326  FillDirManyBody(initialMass, masses, finalState);
327 }
328 
330 FillDirThreeBody(G4double initialMass, const std::vector<G4double>& masses,
331  std::vector<G4LorentzVector>& finalState) {
332  if (GetVerboseLevel()>1)
333  G4cout << " >>> " << GetName() << "::FillDirThreeBody" << G4endl;
334 
335  finalState.resize(3);
336 
337  G4double costh = GenerateCosTheta(kinds[2], modules[2]);
338  finalState[2] = generateWithFixedTheta(costh, modules[2], masses[2]);
339  finalState[2] = toSCM.rotate(finalState[2]); // Align target axis
340 
341  // Generate direction of first particle
342  costh = -0.5 * (modules[2]*modules[2] + modules[0]*modules[0] -
343  modules[1]*modules[1]) / modules[2] / modules[0];
344 
345  if (std::fabs(costh) >= maxCosTheta) { // Bad kinematics; abort generation
346  finalState.clear();
347  return;
348  }
349 
350  // Report success
351  if (GetVerboseLevel()>2) G4cout << " ok for mult 3" << G4endl;
352 
353  // First particle is at fixed angle to recoil system
354  finalState[0] = generateWithFixedTheta(costh, modules[0], masses[0]);
355  finalState[0] = toSCM.rotate(finalState[2], finalState[0]);
356 
357  // Remaining particle is constrained to recoil from entire rest of system
358  finalState[1].set(0.,0.,0.,initialMass);
359  finalState[1] -= finalState[0] + finalState[2];
360 }
361 
363 FillDirManyBody(G4double initialMass, const std::vector<G4double>& masses,
364  std::vector<G4LorentzVector>& finalState) {
365  if (GetVerboseLevel()>1)
366  G4cout << " >>> " << GetName() << "::FillDirManyBody" << G4endl;
367 
368  // Fill all but the last two particles randomly
369  G4double costh = 0.;
370 
371  finalState.resize(multiplicity);
372 
373  for (G4int i=0; i<multiplicity-2; i++) {
374  costh = GenerateCosTheta(kinds[i], modules[i]);
375  finalState[i] = generateWithFixedTheta(costh, modules[i], masses[i]);
376  finalState[i] = toSCM.rotate(finalState[i]); // Align target axis
377  }
378 
379  // Total momentum so far, to compute recoil of last two particles
380  G4LorentzVector psum =
381  std::accumulate(finalState.begin(), finalState.end()-2, G4LorentzVector());
382  G4double pmod = psum.rho();
383 
384  costh = -0.5 * (pmod*pmod +
385  modules[multiplicity-2]*modules[multiplicity-2] -
386  modules[multiplicity-1]*modules[multiplicity-1])
387  / pmod / modules[multiplicity-2];
388 
389  if (GetVerboseLevel() > 2) G4cout << " ct last " << costh << G4endl;
390 
391  if (std::fabs(costh) >= maxCosTheta) { // Bad kinematics; abort generation
392  finalState.clear();
393  return;
394  }
395 
396  // Report success
397  if (GetVerboseLevel()>2) G4cout << " ok for mult " << multiplicity << G4endl;
398 
399  // First particle is at fixed angle to recoil system
400  finalState[multiplicity-2] =
401  generateWithFixedTheta(costh, modules[multiplicity-2],
402  masses[multiplicity-2]);
403  finalState[multiplicity-2] = toSCM.rotate(psum, finalState[multiplicity-2]);
404 
405  // Remaining particle is constrained to recoil from entire rest of system
406  finalState[multiplicity-1].set(0.,0.,0.,initialMass);
407  finalState[multiplicity-1] -= psum + finalState[multiplicity-2];
408 }
409 
410 
411 // Generate polar angle for three- and multi-body systems
412 
414 GenerateCosTheta(G4int ptype, G4double pmod) const {
415  if (GetVerboseLevel() > 2) {
416  G4cout << " >>> " << GetName() << "::GenerateCosTheta " << ptype
417  << " " << pmod << G4endl;
418  }
419 
420  if (multiplicity == 3) { // Use distribution for three-body
421  return angDist->GetCosTheta(bullet_ekin, ptype);
422  }
423 
424  // Throw multi-body distribution
425  G4double p0 = ptype<3 ? 0.36 : 0.25; // Nucleon vs. everything else
426  G4double alf = 1.0 / p0 / (p0 - (pmod+p0)*G4Exp(-pmod / p0));
427 
428  G4double sinth = 2.0;
429 
430  G4int itry1 = -1; /* Loop checking 08.06.2015 MHK */
431  while (std::fabs(sinth) > maxCosTheta && ++itry1 < itry_max) {
432  G4double s1 = pmod * inuclRndm();
433  G4double s2 = alf * oneOverE * p0 * inuclRndm();
434  G4double salf = s1 * alf * G4Exp(-s1 / p0);
435  if (GetVerboseLevel() > 3) {
436  G4cout << " s1 * alf * G4Exp(-s1 / p0) " << salf
437  << " s2 " << s2 << G4endl;
438  }
439 
440  if (salf > s2) sinth = s1 / pmod;
441  }
442 
443  if (GetVerboseLevel() > 3)
444  G4cout << " itry1 " << itry1 << " sinth " << sinth << G4endl;
445 
446  if (itry1 == itry_max) {
447  if (GetVerboseLevel() > 2)
448  G4cout << " high energy angles generation: itry1 " << itry1 << G4endl;
449 
450  sinth = 0.5 * inuclRndm();
451  }
452 
453  // Convert generated sin(theta) to cos(theta) with random sign
454  G4double costh = std::sqrt(1.0 - sinth * sinth);
455  if (inuclRndm() > 0.5) costh = -costh;
456 
457  return costh;
458 }
459 
460 
461 // SPECIAL: Generate N-body phase space using Kopylov algorithm
462 // ==> Code is copied verbatim from G4HadPhaseSpaceKopylov
463 
466  const std::vector<G4double>& masses,
467  std::vector<G4LorentzVector>& finalState) {
468  if (GetVerboseLevel()>2)
469  G4cout << " >>> " << GetName() << "::FillUsingKopylov" << G4endl;
470 
471  finalState.clear();
472 
473  size_t N = masses.size();
474  finalState.resize(N);
475 
476  G4double mtot = std::accumulate(masses.begin(), masses.end(), 0.0);
477  G4double mu = mtot;
478  G4double Mass = initialMass;
479  G4double T = Mass-mtot;
480  G4double recoilMass = 0.0;
481  G4ThreeVector momV, boostV; // Buffers to reduce memory churn
482  G4LorentzVector recoil(0.0,0.0,0.0,Mass);
483 
484  for (size_t k=N-1; k>0; --k) {
485  mu -= masses[k];
486  T *= (k>1) ? BetaKopylov(k) : 0.;
487 
488  recoilMass = mu + T;
489 
490  boostV = recoil.boostVector(); // Previous system's rest frame
491 
492  // Create momentum with a random direction isotropically distributed
493  // FIXME: Should theta distribution should use Bertini fit function?
494  momV.setRThetaPhi(TwoBodyMomentum(Mass,masses[k],recoilMass),
495  UniformTheta(), UniformPhi());
496 
497  finalState[k].setVectM(momV,masses[k]);
498  recoil.setVectM(-momV,recoilMass);
499 
500  finalState[k].boost(boostV);
501  recoil.boost(boostV);
502  Mass = recoilMass;
503  }
504 
505  finalState[0] = recoil;
506 }
507 
509  G4Pow* g4pow = G4Pow::GetInstance();
510 
511  G4int N = 3*K - 5;
512  G4double xN = G4double(N);
513  G4double Fmax = std::sqrt(g4pow->powN(xN/(xN+1.),N)/(xN+1.));
514 
515  G4double F, chi;
516  do { /* Loop checking 08.06.2015 MHK */
517  chi = G4UniformRand();
518  F = std::sqrt(g4pow->powN(chi,N)*(1.-chi));
519  } while ( Fmax*G4UniformRand() > F);
520  return chi;
521 }
virtual const G4String & GetName() const
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
static void setVerboseLevel(G4int vb=0)
static const G4VTwoBodyAngDst * GetDist(G4int is, G4int fs, G4int kw)
Hep3Vector boostVector() const
virtual void GenerateMultiBody(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
const XML_Char * target
Definition: expat.h:268
G4double UniformTheta() const
G4double powN(G4double x, G4int n) const
Definition: G4Pow.cc:128
Definition: G4Pow.hh:56
double x() const
void FillDirThreeBody(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
G4int GetVerboseLevel() const
G4LorentzVector rotate(const G4LorentzVector &mom) const
void setBullet(const G4InuclParticle *bullet)
int G4int
Definition: G4Types.hh:78
double z() const
void SaveKinematics(G4InuclElementaryParticle *bullet, G4InuclElementaryParticle *target)
void setVectM(const Hep3Vector &spatial, double mass)
void FillMagnitudes(G4double initialMass, const std::vector< G4double > &masses)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
void setVerbose(G4int vb=0)
bool G4bool
Definition: G4Types.hh:79
HepLorentzVector & boost(double, double, double)
double rho() const
G4bool satisfyTriangle(const std::vector< G4double > &pmod) const
G4double getKinEnergyInTheTRS() const
void setRThetaPhi(double r, double theta, double phi)
G4double GenerateCosTheta(G4int ptype, G4double pmod) const
virtual const G4String & GetName() const
virtual G4double GetMomentum(G4int ptype, const G4double &ekin) const =0
static void setVerboseLevel(G4int vb=0)
static G4bool usePhaseSpace()
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
virtual void SetVerboseLevel(G4int verbose)
void Configure(G4InuclElementaryParticle *bullet, G4InuclElementaryParticle *target, const std::vector< G4int > &particle_kinds)
const G4String & GetName() const
void set(double x, double y, double z, double t)
virtual void GenerateTwoBody(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
G4LorentzVector generateWithFixedTheta(G4double ct, G4double p, G4double mass=0.)
G4double UniformPhi() const
double y() const
G4double TwoBodyMomentum(G4double M0, G4double M1, G4double M2) const
#define G4endl
Definition: G4ios.hh:61
void FillDirections(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
**D E S C R I P T I O N
Definition: HEPEvtcom.cc:77
void FillUsingKopylov(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
double G4double
Definition: G4Types.hh:76
virtual void SetVerboseLevel(G4int verbose)
void FillDirManyBody(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
virtual G4double GetCosTheta(const G4double &ekin, const G4double &pcm) const =0
static const G4VMultiBodyMomDst * GetDist(G4int is, G4int mult)
G4GLOB_DLL std::ostream G4cerr
void setTarget(const G4InuclParticle *target)
CLHEP::HepLorentzVector G4LorentzVector