Geant4_10
G4INCLClusterDecay.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 // INCL++ intra-nuclear cascade model
27 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
28 // Davide Mancusi, CEA
29 // Alain Boudard, CEA
30 // Sylvie Leray, CEA
31 // Joseph Cugnon, University of Liege
32 //
33 #define INCLXX_IN_GEANT4_MODE 1
34 
35 #include "globals.hh"
36 
44 #include "G4INCLClusterDecay.hh"
45 #include "G4INCLParticleTable.hh"
46 #include "G4INCLKinematicsUtils.hh"
47 #include "G4INCLRandom.hh"
48 #include "G4INCLPhaseSpaceDecay.hh"
49 // #include <cassert>
50 #include <algorithm>
51 
52 namespace G4INCL {
53 
54  namespace ClusterDecay {
55 
56  namespace {
57 
59  void twoBodyDecay(Cluster * const c, ClusterDecayType theDecayMode, ParticleList *decayProducts) {
60  Particle *decayParticle = 0;
61  const ThreeVector mom(0.0, 0.0, 0.0);
62  const ThreeVector pos = c->getPosition();
63 
64  // Create the emitted particle
65  switch(theDecayMode) {
66  case ProtonDecay:
67  decayParticle = new Particle(Proton, mom, pos);
68  break;
69  case NeutronDecay:
70  decayParticle = new Particle(Neutron, mom, pos);
71  break;
72  case AlphaDecay:
73  decayParticle = new Cluster(2,4,false);
74  break;
75  default:
76  INCL_ERROR("Unrecognized cluster-decay mode in two-body decay: " << theDecayMode << std::endl
77  << c->print());
78  return;
79  }
80  decayParticle->makeParticipant();
81  decayParticle->setNumberOfDecays(1);
82  decayParticle->setPosition(c->getPosition());
83  decayParticle->setEmissionTime(c->getEmissionTime());
84  decayParticle->setRealMass();
85 
86  // Save some variables of the mother cluster
87  G4double motherMass = c->getMass();
88  const ThreeVector velocity = -c->boostVector();
89 
90  // Characteristics of the daughter particle
91  const G4int daughterZ = c->getZ() - decayParticle->getZ();
92  const G4int daughterA = c->getA() - decayParticle->getA();
93  const G4double daughterMass = ParticleTable::getRealMass(daughterA,daughterZ);
94 
95  // The mother cluster becomes the daughter
96  c->setZ(daughterZ);
97  c->setA(daughterA);
98  c->setMass(daughterMass);
99  c->setExcitationEnergy(0.);
100 
101  // Decay kinematics in the mother rest frame
102  const G4double decayMass = decayParticle->getMass();
103 // assert(motherMass-daughterMass-decayMass>-1.e-5); // Q-value should be >0
104  G4double pCM = 0.;
105  if(motherMass-daughterMass-decayMass>0.)
106  pCM = KinematicsUtils::momentumInCM(motherMass, daughterMass, decayMass);
107  const ThreeVector momentum = Random::normVector(pCM);
108  c->setMomentum(momentum);
109  c->adjustEnergyFromMomentum();
110  decayParticle->setMomentum(-momentum);
111  decayParticle->adjustEnergyFromMomentum();
112 
113  // Boost to the lab frame
114  decayParticle->boost(velocity);
115  c->boost(velocity);
116 
117  // Add the decay particle to the list of decay products
118  decayProducts->push_back(decayParticle);
119  }
120 
122  void threeBodyDecay(Cluster * const c, ClusterDecayType theDecayMode, ParticleList *decayProducts) {
123  Particle *decayParticle1 = 0;
124  Particle *decayParticle2 = 0;
125  const ThreeVector mom(0.0, 0.0, 0.0);
126  const ThreeVector pos = c->getPosition();
127 
128  // Create the emitted particles
129  switch(theDecayMode) {
130  case TwoProtonDecay:
131  decayParticle1 = new Particle(Proton, mom, pos);
132  decayParticle2 = new Particle(Proton, mom, pos);
133  break;
134  case TwoNeutronDecay:
135  decayParticle1 = new Particle(Neutron, mom, pos);
136  decayParticle2 = new Particle(Neutron, mom, pos);
137  break;
138  default:
139  INCL_ERROR("Unrecognized cluster-decay mode in three-body decay: " << theDecayMode << std::endl
140  << c->print());
141  return;
142  }
143  decayParticle1->makeParticipant();
144  decayParticle2->makeParticipant();
145  decayParticle1->setNumberOfDecays(1);
146  decayParticle2->setNumberOfDecays(1);
147  decayParticle1->setRealMass();
148  decayParticle2->setRealMass();
149 
150  // Save some variables of the mother cluster
151  const G4double motherMass = c->getMass();
152  const ThreeVector velocity = -c->boostVector();
153 
154  // Masses and charges of the daughter particle and of the decay products
155  const G4int decayZ1 = decayParticle1->getZ();
156  const G4int decayA1 = decayParticle1->getA();
157  const G4int decayZ2 = decayParticle2->getZ();
158  const G4int decayA2 = decayParticle2->getA();
159  const G4int decayZ = decayZ1 + decayZ2;
160  const G4int decayA = decayA1 + decayA2;
161  const G4int daughterZ = c->getZ() - decayZ;
162  const G4int daughterA = c->getA() - decayA;
163  const G4double decayMass1 = decayParticle1->getMass();
164  const G4double decayMass2 = decayParticle2->getMass();
165  const G4double daughterMass = ParticleTable::getRealMass(daughterA,daughterZ);
166 
167  // Q-values
168  G4double qValue = motherMass - daughterMass - decayMass1 - decayMass2;
169 // assert(qValue > -1e-5); // Q-value should be >0
170  if(qValue<0.)
171  qValue=0.;
172  const G4double qValueB = qValue * Random::shoot();
173 
174  // The decay particles behave as if they had more mass until the second
175  // decay
176  const G4double decayMass = decayMass1 + decayMass2 + qValueB;
177 
178  /* Stage A: mother --> daughter + (decay1+decay2) */
179 
180  // The mother cluster becomes the daughter
181  c->setZ(daughterZ);
182  c->setA(daughterA);
183  c->setMass(daughterMass);
184  c->setExcitationEnergy(0.);
185 
186  // Decay kinematics in the mother rest frame
187  const G4double pCMA = KinematicsUtils::momentumInCM(motherMass, daughterMass, decayMass);
188  const ThreeVector momentumA = Random::normVector(pCMA);
189  c->setMomentum(momentumA);
190  c->adjustEnergyFromMomentum();
191  const ThreeVector decayBoostVector = momentumA/std::sqrt(decayMass*decayMass + momentumA.mag2());
192 
193  /* Stage B: (decay1+decay2) --> decay1 + decay2 */
194 
195  // Decay kinematics in the (decay1+decay2) rest frame
196  const G4double pCMB = KinematicsUtils::momentumInCM(decayMass, decayMass1, decayMass2);
197  const ThreeVector momentumB = Random::normVector(pCMB);
198  decayParticle1->setMomentum(momentumB);
199  decayParticle2->setMomentum(-momentumB);
200  decayParticle1->adjustEnergyFromMomentum();
201  decayParticle2->adjustEnergyFromMomentum();
202 
203  // Boost decay1 and decay2 to the Stage-A decay frame
204  decayParticle1->boost(decayBoostVector);
205  decayParticle2->boost(decayBoostVector);
206 
207  // Boost all particles to the lab frame
208  decayParticle1->boost(velocity);
209  decayParticle2->boost(velocity);
210  c->boost(velocity);
211 
212  // Add the decay particles to the list of decay products
213  decayProducts->push_back(decayParticle1);
214  decayProducts->push_back(decayParticle2);
215  }
216 
217 #ifndef INCLXX_IN_GEANT4_MODE
218 
230  void phaseSpaceDecayLegacy(Cluster * const c, ClusterDecayType theDecayMode, ParticleList *decayProducts) {
231  const G4int theA = c->getA();
232  const G4int theZ = c->getZ();
233  const ThreeVector mom(0.0, 0.0, 0.0);
234  const ThreeVector pos = c->getPosition();
235 
236  G4int theZStep;
237  ParticleType theEjectileType;
238  switch(theDecayMode) {
239  case ProtonUnbound:
240  theZStep = 1;
241  theEjectileType = Proton;
242  break;
243  case NeutronUnbound:
244  theZStep = 0;
245  theEjectileType = Neutron;
246  break;
247  default:
248  INCL_ERROR("Unrecognized cluster-decay mode in phase-space decay: " << theDecayMode << std::endl
249  << c->print());
250  return;
251  }
252 
253  // Find the daughter cluster (first cluster which is not
254  // proton/neutron-unbound, in the sense of the table)
255  G4int finalDaughterZ, finalDaughterA;
257  finalDaughterZ=theZ;
258  finalDaughterA=theA;
259  while(clusterDecayMode[finalDaughterZ][finalDaughterA]==theDecayMode) {
260  finalDaughterA--;
261  finalDaughterZ -= theZStep;
262  }
263  } else {
264  finalDaughterA = 1;
265  if(theDecayMode==ProtonUnbound)
266  finalDaughterZ = 1;
267  else
268  finalDaughterZ = 0;
269  }
270 // assert(finalDaughterZ<=theZ && finalDaughterA<theA && finalDaughterA>0 && finalDaughterZ>=0);
271  const G4double finalDaughterMass = ParticleTable::getRealMass(finalDaughterA, finalDaughterZ);
272 
273  // Compute the available decay energy
274  const G4int nSplits = theA-finalDaughterA;
275  const G4double ejectileMass = ParticleTable::getRealMass(1, theZStep);
276  // c->getMass() can possibly contain some excitation energy, too
277  G4double availableEnergy = c->getMass() - finalDaughterMass - nSplits*ejectileMass;
278 // assert(availableEnergy>-1.e-5);
279  if(availableEnergy<0.)
280  availableEnergy=0.;
281 
282  // Compute an estimate of the maximum event weight
283  G4double maximumWeight = 1.;
284  G4double eMax = finalDaughterMass + availableEnergy;
285  G4double eMin = finalDaughterMass - ejectileMass;
286  for(G4int iSplit=0; iSplit<nSplits; ++iSplit) {
287  eMax += ejectileMass;
288  eMin += ejectileMass;
289  maximumWeight *= KinematicsUtils::momentumInCM(eMax, eMin, ejectileMass);
290  }
291 
292  // Sample decays until the weight cutoff is satisfied
294  std::vector<G4double> theCMMomenta;
295  std::vector<G4double> invariantMasses;
296  G4int nTries=0;
297  /* Maximum number of trials dependent on nSplits. 50 trials seems to be
298  * sufficient for small nSplits. For nSplits>=5, maximumWeight is a gross
299  * overestimate of the actual maximum weight, leading to unreasonably high
300  * rejection rates. For these cases, we set nSplits=1000, although the sane
301  * thing to do would be to improve the importance sampling (maybe by
302  * parametrising maximumWeight?).
303  */
304  G4int maxTries;
305  if(nSplits<5)
306  maxTries=50;
307  else
308  maxTries=1000;
309  do {
310  if(nTries++>maxTries) {
311  INCL_WARN("Phase-space decay exceeded the maximum number of rejections (" << maxTries
312  << "). Z=" << theZ << ", A=" << theA << ", E*=" << c->getExcitationEnergy()
313  << ", availableEnergy=" << availableEnergy
314  << ", nSplits=" << nSplits
315  << std::endl);
316  break;
317  }
318 
319  // Construct a sorted vector of random numbers
320  std::vector<G4double> randomNumbers;
321  for(G4int iSplit=0; iSplit<nSplits-1; ++iSplit)
322  randomNumbers.push_back(Random::shoot0());
323  std::sort(randomNumbers.begin(), randomNumbers.end());
324 
325  // Divide the available decay energy in the right number of steps
326  invariantMasses.clear();
327  invariantMasses.push_back(finalDaughterMass);
328  for(G4int iSplit=0; iSplit<nSplits-1; ++iSplit)
329  invariantMasses.push_back(finalDaughterMass + (iSplit+1)*ejectileMass + randomNumbers.at(iSplit)*availableEnergy);
330  invariantMasses.push_back(c->getMass());
331 
332  weight = 1.;
333  theCMMomenta.clear();
334  for(G4int iSplit=0; iSplit<nSplits; ++iSplit) {
335  G4double motherMass = invariantMasses.at(nSplits-iSplit);
336  const G4double daughterMass = invariantMasses.at(nSplits-iSplit-1);
337 // assert(motherMass-daughterMass-ejectileMass>-1.e-5);
338  G4double pCM = 0.;
339  if(motherMass-daughterMass-ejectileMass>0.)
340  pCM = KinematicsUtils::momentumInCM(motherMass, daughterMass, ejectileMass);
341  theCMMomenta.push_back(pCM);
342  weight *= pCM;
343  }
344  } while(maximumWeight*Random::shoot()>weight);
345 
346  for(G4int iSplit=0; iSplit<nSplits; ++iSplit) {
347  ThreeVector const velocity = -c->boostVector();
348 
349 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEANT4_MODE)
350  const G4double motherMass = c->getMass();
351 #endif
352  c->setA(c->getA() - 1);
353  c->setZ(c->getZ() - theZStep);
354  c->setMass(invariantMasses.at(nSplits-iSplit-1));
355 
356  Particle *ejectile = new Particle(theEjectileType, mom, pos);
357  ejectile->setRealMass();
358 
359 // assert(motherMass-c->getMass()-ejectileMass>-1.e-5);
360  ThreeVector momentum;
361  momentum = Random::normVector(theCMMomenta.at(iSplit));
362  c->setMomentum(momentum);
363  c->adjustEnergyFromMomentum();
364  ejectile->setMomentum(-momentum);
365  ejectile->adjustEnergyFromMomentum();
366 
367  // Boost to the lab frame
368  c->boost(velocity);
369  ejectile->boost(velocity);
370 
371  // Add the decay particle to the list of decay products
372  decayProducts->push_back(ejectile);
373  }
374 // assert(std::abs(c->getRealMass()-c->getMass())<1.e-3);
375  c->setExcitationEnergy(0.);
376  }
377 #endif // INCLXX_IN_GEANT4_MODE
378 
384  void phaseSpaceDecay(Cluster * const c, ClusterDecayType theDecayMode, ParticleList *decayProducts) {
385  const G4int theA = c->getA();
386  const G4int theZ = c->getZ();
387  const ThreeVector mom(0.0, 0.0, 0.0);
388  const ThreeVector pos = c->getPosition();
389 
390  G4int theZStep;
391  ParticleType theEjectileType;
392  switch(theDecayMode) {
393  case ProtonUnbound:
394  theZStep = 1;
395  theEjectileType = Proton;
396  break;
397  case NeutronUnbound:
398  theZStep = 0;
399  theEjectileType = Neutron;
400  break;
401  default:
402  INCL_ERROR("Unrecognized cluster-decay mode in phase-space decay: " << theDecayMode << std::endl
403  << c->print());
404  return;
405  }
406 
407  // Find the daughter cluster (first cluster which is not
408  // proton/neutron-unbound, in the sense of the table)
409  G4int finalDaughterZ, finalDaughterA;
411  finalDaughterZ=theZ;
412  finalDaughterA=theA;
413  while(finalDaughterA>0 && clusterDecayMode[finalDaughterZ][finalDaughterA]==theDecayMode) {
414  finalDaughterA--;
415  finalDaughterZ -= theZStep;
416  }
417  } else {
418  finalDaughterA = 1;
419  if(theDecayMode==ProtonUnbound)
420  finalDaughterZ = 1;
421  else
422  finalDaughterZ = 0;
423  }
424 // assert(finalDaughterZ<=theZ && finalDaughterA<theA && finalDaughterA>0 && finalDaughterZ>=0);
425 
426  // Compute the available decay energy
427  const G4int nSplits = theA-finalDaughterA;
428  // c->getMass() can possibly contain some excitation energy, too
429  const G4double availableEnergy = c->getMass();
430 
431  // Save the boost vector for the cluster
432  const ThreeVector boostVector = - c->boostVector();
433 
434  // Create a list of decay products
435  ParticleList products;
436  c->setA(finalDaughterA);
437  c->setZ(finalDaughterZ);
438  c->setRealMass();
439  c->setMomentum(ThreeVector());
440  c->adjustEnergyFromMomentum();
441  products.push_back(c);
442  for(G4int i=0; i<nSplits; ++i) {
443  Particle *ejectile = new Particle(theEjectileType, mom, pos);
444  ejectile->setRealMass();
445  products.push_back(ejectile);
446  }
447 
448  G4INCL::PhaseSpaceDecay::decay(availableEnergy, boostVector, products);
449 
450  // Copy decay products in the output list (but skip the residue)
451  ParticleList::iterator productsIter = products.begin();
452  std::advance(productsIter, 1);
453  decayProducts->insert(decayProducts->end(), productsIter, products.end());
454 
455  c->setExcitationEnergy(0.);
456  }
457 
463  void recursiveDecay(Cluster * const c, ParticleList *decayProducts) {
464  const G4int Z = c->getZ();
465  const G4int A = c->getA();
466 // assert(c->getExcitationEnergy()>-1.e-5);
467  if(c->getExcitationEnergy()<0.)
468  c->setExcitationEnergy(0.);
469 
471  ClusterDecayType theDecayMode = clusterDecayMode[Z][A];
472 
473  switch(theDecayMode) {
474  default:
475  INCL_ERROR("Unrecognized cluster-decay mode: " << theDecayMode << std::endl
476  << c->print());
477  return;
478  break;
479  case StableCluster:
480  // For stable clusters, just return
481  return;
482  break;
483  case ProtonDecay:
484  case NeutronDecay:
485  case AlphaDecay:
486  // Two-body decays
487  twoBodyDecay(c, theDecayMode, decayProducts);
488  break;
489  case TwoProtonDecay:
490  case TwoNeutronDecay:
491  // Three-body decays
492  threeBodyDecay(c, theDecayMode, decayProducts);
493  break;
494  case ProtonUnbound:
495  case NeutronUnbound:
496  // Phase-space decays
497  phaseSpaceDecay(c, theDecayMode, decayProducts);
498  break;
499  }
500 
501  // Calls itself recursively in case the produced remnant is still unstable.
502  // Sneaky, isn't it.
503  recursiveDecay(c,decayProducts);
504 
505  } else {
506  // The cluster is too large for our decay-mode table. Decompose it only
507  // if Z==0 || Z==A.
508  INCL_DEBUG("Cluster is outside the decay-mode table." << c->print() << std::endl);
509  if(Z==A) {
510  INCL_DEBUG("Z==A, will decompose it in free protons." << std::endl);
511  phaseSpaceDecay(c, ProtonUnbound, decayProducts);
512  } else if(Z==0) {
513  INCL_DEBUG("Z==0, will decompose it in free neutrons." << std::endl);
514  phaseSpaceDecay(c, NeutronUnbound, decayProducts);
515  }
516  }
517  }
518 
519  } // namespace
520 
521  G4bool isStable(Cluster const * const c) {
522  const G4int Z = c->getZ();
523  const G4int A = c->getA();
524  return (clusterDecayMode[Z][A]==StableCluster);
525  }
526 
538  {
539  /* A = 0 1 2 3 4 5 6 7 8 9 10 11 12 */
549  };
550 
552  ParticleList decayProducts;
553  recursiveDecay(c, &decayProducts);
554 
555  // Correctly update the particle type
556  if(c->getA()==1) {
557 // assert(c->getZ()==1 || c->getZ()==0);
558  if(c->getZ()==1)
559  c->setType(Proton);
560  else
561  c->setType(Neutron);
562  c->setRealMass();
563  }
564 
565  return decayProducts;
566  }
567 
568  } // namespace ClusterDecay
569 } // namespace G4INCL
570 
G4int getA() const
Returns the baryon number.
Static class for carrying out cluster decays.
#define INCL_ERROR(x)
G4ThreadLocal ClusterDecayType clusterDecayMode[ParticleTable::clusterTableZSize][ParticleTable::clusterTableASize]
Table for cluster decays.
#define INCL_WARN(x)
G4bool isStable(Cluster const *const c)
True if the cluster is stable.
subroutine sort(A, N)
Definition: dpm25nuc7.f:4670
G4double momentumInCM(Particle const *const p1, Particle const *const p2)
gives the momentum in the CM frame of two particles.
#define G4ThreadLocal
Definition: tls.hh:52
double weight
Definition: plottest35.C:25
G4double shoot0()
Definition: G4INCLRandom.cc:78
int G4int
Definition: G4Types.hh:78
ThreeVector normVector(G4double norm=1.)
G4double getRealMass(const G4INCL::ParticleType t)
Get particle mass (in MeV/c^2)
UnorderedVector< Particle * > ParticleList
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
void decay(G4double initialMass, const ThreeVector &theBoostVector, ParticleList &particles)
Generate decay momenta according to a uniform phase-space model.
Float_t Z
Definition: plot.C:39
bool G4bool
Definition: G4Types.hh:79
G4int getZ() const
Returns the charge number.
void setRealMass()
Set the mass of the Particle to its real mass.
void setType(ParticleType t)
G4double shoot()
Definition: G4INCLRandom.cc:74
double G4double
Definition: G4Types.hh:76
#define INCL_DEBUG(x)
tuple c
Definition: test.py:13
G4double getZ() const