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