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