| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 H
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
38 /** \file G4INCLClusterDecay.cc
39 * \brief Static class for carrying out cluste
40 *
41 * \date 6th July 2011
42 * \author Davide Mancusi
43 */
44
45 #include "G4INCLClusterDecay.hh"
46 #include "G4INCLParticleTable.hh"
47 #include "G4INCLKinematicsUtils.hh"
48 #include "G4INCLRandom.hh"
49 #include "G4INCLPhaseSpaceGenerator.hh"
50 // #include <cassert>
51 #include <algorithm>
52
53 namespace G4INCL {
54
55 namespace ClusterDecay {
56
57 namespace {
58
59 /// \brief Carries out two-body decays
60 void twoBodyDecay(Cluster * const c, Clu
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(Proto
69 break;
70 case NeutronDecay:
71 decayParticle = new Particle(Neutr
72 break;
73 case AlphaDecay:
74 decayParticle = new Cluster(2,4,0,
75 break;
76 case LambdaDecay:
77 decayParticle = new Particle(Lambd
78 break;
79 default:
80 INCL_ERROR("Unrecognized cluster-d
81 << c->print());
82 return;
83 }
84 decayParticle->makeParticipant();
85 decayParticle->setNumberOfDecays(1);
86 decayParticle->setPosition(c->getPosit
87 decayParticle->setEmissionTime(c->getE
88 decayParticle->setRealMass();
89
90 // Save some variables of the mother c
91 #ifdef INCLXX_IN_GEANT4_MODE
92 if ((c->getZ() == 1) && (c->getA() =
93 c->setMass(2053.952);
94 if (c->getEnergy() < 2053.952)
95 c->setMomentum(c->getMomentu
96 else
97 c->setMomentum(c->getMomentu
98 }
99 #endif
100 G4double motherMass = c->getMass();
101 const ThreeVector velocity = -c->boost
102
103 // Characteristics of the daughter par
104 const G4int daughterZ = c->getZ() - de
105 const G4int daughterA = c->getA() - de
106 const G4int daughterS = c->getS() - de
107 const G4double daughterMass = Particle
108
109 // The mother cluster becomes the daug
110 c->setZ(daughterZ);
111 c->setA(daughterA);
112 c->setS(daughterS);
113 c->setMass(daughterMass);
114 c->setExcitationEnergy(0.);
115
116 // Decay kinematics in the mother rest
117 const G4double decayMass = decayPartic
118 // assert(motherMass-daughterMass-decayMass>-1
119 G4double pCM = 0.;
120 if(motherMass-daughterMass-decayMass>0
121 pCM = KinematicsUtils::momentumInCM(
122 const ThreeVector momentum = Random::n
123 c->setMomentum(momentum);
124 c->adjustEnergyFromMomentum();
125 decayParticle->setMomentum(-momentum);
126 decayParticle->adjustEnergyFromMomentu
127
128 // Boost to the lab frame
129 decayParticle->boost(velocity);
130 c->boost(velocity);
131
132 // Add the decay particle to the list
133 decayProducts->push_back(decayParticle
134 }
135
136 /// \brief Carries out three-body decays
137 void threeBodyDecay(Cluster * const c, C
138 Particle *decayParticle1 = 0;
139 Particle *decayParticle2 = 0;
140 const ThreeVector mom(0.0, 0.0, 0.0);
141 const ThreeVector pos = c->getPosition
142
143 // Create the emitted particles
144 switch(theDecayMode) {
145 case TwoProtonDecay:
146 decayParticle1 = new Particle(Prot
147 decayParticle2 = new Particle(Prot
148 break;
149 case TwoNeutronDecay:
150 decayParticle1 = new Particle(Neut
151 decayParticle2 = new Particle(Neut
152 break;
153 default:
154 INCL_ERROR("Unrecognized cluster-d
155 << c->print());
156 return;
157 }
158 decayParticle1->makeParticipant();
159 decayParticle2->makeParticipant();
160 decayParticle1->setNumberOfDecays(1);
161 decayParticle2->setNumberOfDecays(1);
162 decayParticle1->setRealMass();
163 decayParticle2->setRealMass();
164
165 // Save some variables of the mother c
166 const G4double motherMass = c->getMass
167 const ThreeVector velocity = -c->boost
168
169 // Masses and charges of the daughter
170 const G4int decayZ1 = decayParticle1->
171 const G4int decayA1 = decayParticle1->
172 const G4int decayS1 = decayParticle1->
173 const G4int decayZ2 = decayParticle2->
174 const G4int decayA2 = decayParticle2->
175 const G4int decayS2 = decayParticle2->
176 const G4int decayZ = decayZ1 + decayZ2
177 const G4int decayA = decayA1 + decayA2
178 const G4int decayS = decayS1 + decayS2
179 const G4int daughterZ = c->getZ() - de
180 const G4int daughterA = c->getA() - de
181 const G4int daughterS = c->getS() - de
182 const G4double decayMass1 = decayParti
183 const G4double decayMass2 = decayParti
184 const G4double daughterMass = Particle
185
186 // Q-values
187 G4double qValue = motherMass - daughte
188 // assert(qValue > -1e-5); // Q-value should b
189 if(qValue<0.)
190 qValue=0.;
191 const G4double qValueB = qValue * Rand
192
193 // The decay particles behave as if th
194 // decay
195 const G4double decayMass = decayMass1
196
197 /* Stage A: mother --> daughter + (dec
198
199 // The mother cluster becomes the daug
200 c->setZ(daughterZ);
201 c->setA(daughterA);
202 c->setS(daughterS);
203 c->setMass(daughterMass);
204 c->setExcitationEnergy(0.);
205
206 // Decay kinematics in the mother rest
207 const G4double pCMA = KinematicsUtils:
208 const ThreeVector momentumA = Random::
209 c->setMomentum(momentumA);
210 c->adjustEnergyFromMomentum();
211 const ThreeVector decayBoostVector = m
212
213 /* Stage B: (decay1+decay2) --> decay1
214
215 // Decay kinematics in the (decay1+dec
216 const G4double pCMB = KinematicsUtils:
217 const ThreeVector momentumB = Random::
218 decayParticle1->setMomentum(momentumB)
219 decayParticle2->setMomentum(-momentumB
220 decayParticle1->adjustEnergyFromMoment
221 decayParticle2->adjustEnergyFromMoment
222
223 // Boost decay1 and decay2 to the Stag
224 decayParticle1->boost(decayBoostVector
225 decayParticle2->boost(decayBoostVector
226
227 // Boost all particles to the lab fram
228 decayParticle1->boost(velocity);
229 decayParticle2->boost(velocity);
230 c->boost(velocity);
231
232 // Add the decay particles to the list
233 decayProducts->push_back(decayParticle
234 decayProducts->push_back(decayParticle
235 }
236
237 #ifdef INCL_DO_NOT_COMPILE
238 /** \brief Disassembles unbound nuclei u
239 *
240 * The decay products are assumed to uni
241 * (the "phase-space" naming is a long-s
242 * The generation of the final state is
243 * Raubold-Lynch method. Parts of our im
244 * ROOT's TGenPhaseSpace class, which in
245 * historical GENBOD routine [CERN repor
246 * implementation is documented at the f
247 *
248 * http://root.cern.ch/root/html/TGenPha
249 */
250 void phaseSpaceDecayLegacy(Cluster * con
251 const G4int theA = c->getA();
252 const G4int theZ = c->getZ();
253 // assert(c->getS() == 0);
254 const ThreeVector mom(0.0, 0.0, 0.0);
255 const ThreeVector pos = c->getPosition
256
257 G4int theZStep;
258 ParticleType theEjectileType;
259 switch(theDecayMode) {
260 case ProtonUnbound:
261 theZStep = 1;
262 theEjectileType = Proton;
263 break;
264 case NeutronUnbound:
265 theZStep = 0;
266 theEjectileType = Neutron;
267 break;
268 default:
269 INCL_ERROR("Unrecognized cluster-d
270 << c->print());
271 return;
272 }
273
274 // Find the daughter cluster (first cl
275 // proton/neutron-unbound, in the sens
276 G4int finalDaughterZ, finalDaughterA;
277 if(theZ<ParticleTable::clusterTableZSi
278 finalDaughterZ=theZ;
279 finalDaughterA=theA;
280 while(clusterDecayMode[0][finalDaugh
281 finalDaughterA--;
282 finalDaughterZ -= theZStep;
283 }
284 } else {
285 finalDaughterA = 1;
286 if(theDecayMode==ProtonUnbound)
287 finalDaughterZ = 1;
288 else
289 finalDaughterZ = 0;
290 }
291 // assert(finalDaughterZ<=theZ && finalDaughte
292 const G4double finalDaughterMass = Par
293
294 // Compute the available decay energy
295 const G4int nSplits = theA-finalDaught
296 const G4double ejectileMass = Particle
297 // c->getMass() can possibly contain s
298 G4double availableEnergy = c->getMass(
299 // assert(availableEnergy>-1.e-5);
300 if(availableEnergy<0.)
301 availableEnergy=0.;
302
303 // Compute an estimate of the maximum
304 G4double maximumWeight = 1.;
305 G4double eMax = finalDaughterMass + av
306 G4double eMin = finalDaughterMass - ej
307 for(G4int iSplit=0; iSplit<nSplits; ++
308 eMax += ejectileMass;
309 eMin += ejectileMass;
310 maximumWeight *= KinematicsUtils::mo
311 }
312
313 // Sample decays until the weight cuto
314 G4double weight;
315 std::vector<G4double> theCMMomenta;
316 std::vector<G4double> invariantMasses;
317 G4int nTries=0;
318 /* Maximum number of trials dependent
319 * sufficient for small nSplits. For n
320 * overestimate of the actual maximum
321 * rejection rates. For these cases, w
322 * thing to do would be to improve the
323 * parametrising maximumWeight?).
324 */
325 G4int maxTries;
326 if(nSplits<5)
327 maxTries=50;
328 else
329 maxTries=1000;
330 do {
331 if(nTries++>maxTries) {
332 INCL_WARN("Phase-space decay excee
333 << "). Z=" << theZ << ", A="
334 << ", availableEnergy=" << av
335 << ", nSplits=" << nSplits
336 << '\n');
337 break;
338 }
339
340 // Construct a sorted vector of rand
341 std::vector<G4double> randomNumbers;
342 for(G4int iSplit=0; iSplit<nSplits-1
343 randomNumbers.push_back(Random::sh
344 std::sort(randomNumbers.begin(), ran
345
346 // Divide the available decay energy
347 invariantMasses.clear();
348 invariantMasses.push_back(finalDaugh
349 for(G4int iSplit=0; iSplit<nSplits-1
350 invariantMasses.push_back(finalDau
351 invariantMasses.push_back(c->getMass
352
353 weight = 1.;
354 theCMMomenta.clear();
355 for(G4int iSplit=0; iSplit<nSplits;
356 G4double motherMass = invariantMas
357 const G4double daughterMass = inva
358 // assert(motherMass-daughterMass-ejectileMass
359 G4double pCM = 0.;
360 if(motherMass-daughterMass-ejectil
361 pCM = KinematicsUtils::momentumI
362 theCMMomenta.push_back(pCM);
363 weight *= pCM;
364 }
365 } while(maximumWeight*Random::shoot()>
366
367 for(G4int iSplit=0; iSplit<nSplits; ++
368 ThreeVector const velocity = -c->boo
369
370 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEA
371 const G4double motherMass = c->getMa
372 #endif
373 c->setA(c->getA() - 1);
374 c->setZ(c->getZ() - theZStep);
375 c->setMass(invariantMasses.at(nSplit
376
377 Particle *ejectile = new Particle(th
378 ejectile->setRealMass();
379
380 // assert(motherMass-c->getMass()-ejectileMass
381 ThreeVector momentum;
382 momentum = Random::normVector(theCMM
383 c->setMomentum(momentum);
384 c->adjustEnergyFromMomentum();
385 ejectile->setMomentum(-momentum);
386 ejectile->adjustEnergyFromMomentum()
387
388 // Boost to the lab frame
389 c->boost(velocity);
390 ejectile->boost(velocity);
391
392 // Add the decay particle to the lis
393 decayProducts->push_back(ejectile);
394 }
395 // assert(std::abs(c->getRealMass()-c->getMass
396 c->setExcitationEnergy(0.);
397 }
398 #endif // INCL_DO_NOT_COMPILE
399
400 /** \brief Disassembles unbound nuclei u
401 *
402 * This implementation uses the Kopylov
403 * PhaseSpaceGenerator.
404 */
405 void phaseSpaceDecay(Cluster * const c,
406 const G4int theA = c->getA();
407 const G4int theZ = c->getZ();
408 const G4int theL = (-1)*(c->getS());
409 const ThreeVector mom(0.0, 0.0, 0.0);
410 const ThreeVector pos = c->getPosition
411
412 G4int theZStep;
413
414 ParticleType theEjectileType;
415 switch(theDecayMode) {
416 case ProtonUnbound:
417 theZStep = 1;
418 theEjectileType = Proton;
419 break;
420 case NeutronUnbound:
421 theZStep = 0;
422 theEjectileType = Neutron;
423 break;
424 case LambdaUnbound: // Will always c
425 theZStep = -99;
426 if(theZ==0) theEjectileType = Neut
427 else theEjectileType = Proton;
428 break;
429 default:
430 INCL_ERROR("Unrecognized cluster-d
431 << c->print());
432 return;
433 }
434
435 // Find the daughter cluster (first cl
436 // proton/neutron-unbound, in the sens
437 G4int finalDaughterZ, finalDaughterA,
438 if(theZ<ParticleTable::clusterTableZSi
439 finalDaughterZ=theZ;
440 finalDaughterA=theA;
441 finalDaughterL=theL;
442 while(finalDaughterA>0 && clusterDec
443 finalDaughterA--;
444 finalDaughterZ -= theZStep;
445 }
446 } else {
447 finalDaughterA = 1;
448 if(theDecayMode==ProtonUnbound){
449 finalDaughterZ = 1;
450 finalDaughterL = 0;
451 }
452 else if(theDecayMode==NeutronUnbound
453 finalDaughterZ = 0;
454 finalDaughterL = 0;
455 }
456 else {
457 finalDaughterZ = 0;
458 finalDaughterL = 1;
459 }
460 }
461 // assert(finalDaughterZ<=theZ && finalDaughte
462
463 // Compute the available decay energy
464 const G4int nLambda = theL-finalDaught
465 const G4int nSplits = theA-finalDaught
466 // c->getMass() can possibly contain s
467 const G4double availableEnergy = c->ge
468
469 // Save the boost vector for the clust
470 const ThreeVector boostVector = - c->b
471
472 // Create a list of decay products
473 ParticleList products;
474 c->setA(finalDaughterA);
475 c->setZ(finalDaughterZ);
476 c->setS((-1)*finalDaughterL);
477 c->setRealMass();
478 c->setMomentum(ThreeVector());
479 c->adjustEnergyFromMomentum();
480 products.push_back(c);
481
482 for(G4int j=0; j<nLambda; ++j) {
483 Particle *ejectile = new Particle(Lambda
484 ejectile->setRealMass();
485 products.push_back(ejectile);
486 }
487 for(G4int i=0; i<nSplits; ++i) {
488 Particle *ejectile = new Particle(th
489 ejectile->setRealMass();
490 products.push_back(ejectile);
491 }
492
493 PhaseSpaceGenerator::generate(availabl
494 products.boost(boostVector);
495
496 // Copy decay products in the output l
497 ParticleList::iterator productsIter =
498 std::advance(productsIter, 1);
499 decayProducts->insert(decayProducts->e
500
501 c->setExcitationEnergy(0.);
502 }
503
504 /** \brief Recursively decay clusters
505 *
506 * \param c cluster that should decay
507 * \param decayProducts decay products a
508 */
509 void recursiveDecay(Cluster * const c, P
510 const G4int Z = c->getZ();
511 const G4int A = c->getA();
512 const G4int S = c->getS();
513 // assert(c->getExcitationEnergy()>-1.e-5);
514 if(c->getExcitationEnergy()<0.)
515 c->setExcitationEnergy(0.);
516
517 if(Z<ParticleTable::clusterTableZSize
518 ClusterDecayType theDecayMode = clus
519
520 switch(theDecayMode) {
521 default:
522 INCL_ERROR("Unrecognized cluster
523 << c->print());
524 return;
525 break;
526 case StableCluster:
527 // For stable clusters, just ret
528 return;
529 break;
530 case ProtonDecay:
531 case NeutronDecay:
532 case LambdaDecay:
533 case AlphaDecay:
534 // Two-body decays
535 twoBodyDecay(c, theDecayMode, de
536 break;
537 case TwoProtonDecay:
538 case TwoNeutronDecay:
539 // Three-body decays
540 threeBodyDecay(c, theDecayMode,
541 break;
542 case ProtonUnbound:
543 case NeutronUnbound:
544 case LambdaUnbound:
545 // Phase-space decays
546 phaseSpaceDecay(c, theDecayMode,
547 break;
548 }
549
550 // Calls itself recursively in case
551 // Sneaky, isn't it.
552 recursiveDecay(c,decayProducts);
553
554 } else {
555 // The cluster is too large for our
556 // if Z==0 || Z==A.
557 INCL_DEBUG("Cluster is outside the d
558 if(Z==A) {
559 INCL_DEBUG("Z==A, will decompose i
560 phaseSpaceDecay(c, ProtonUnbound,
561 } else if(Z==0) {
562 INCL_DEBUG("Z==0, will decompose i
563 phaseSpaceDecay(c, NeutronUnbound,
564 }
565 }
566 }
567
568 } // namespace
569
570 G4bool isStable(Cluster const * const c) {
571 const G4int Z = c->getZ();
572 const G4int A = c->getA();
573 const G4int L = ((-1)*c->getS());
574 return (clusterDecayMode[L][Z][A]==Stabl
575 }
576
577 /** \brief Table for cluster decays
578 *
579 * Definition of "Stable": halflife > 1 ms
580 *
581 * These table includes decay data for clu
582 * not produce. It can't hurt.
583 *
584 * Unphysical nuclides (A<Z) are marked as
585 * produced by INCL. If you find them in t
586 */
587 G4ThreadLocal ClusterDecayType clusterDeca
588 {{/* S = 0 */
589 /* A = 0
590 /* Z = 0 */ {StableCluster, StableCl
591 /* Z = 1 */ {StableCluster, StableCl
592 /* Z = 2 */ {StableCluster, StableCl
593 /* Z = 3 */ {StableCluster, StableCl
594 /* Z = 4 */ {StableCluster, StableCl
595 /* Z = 5 */ {StableCluster, StableCl
596 /* Z = 6 */ {StableCluster, StableCl
597 /* Z = 7 */ {StableCluster, StableCl
598 /* Z = 8 */ {StableCluster, StableCl
599 },
600 { /* S = -1 */
601 /* A = 0
602 /* Z = 0 */ {StableCluster, StableCl
603 /* Z = 1 */ {StableCluster, StableCl
604 /* Z = 2 */ {StableCluster, StableCl
605 /* Z = 3 */ {StableCluster, StableCl
606 /* Z = 4 */ {StableCluster, StableCl
607 /* Z = 5 */ {StableCluster, StableCl
608 /* Z = 6 */ {StableCluster, StableCl
609 /* Z = 7 */ {StableCluster, StableCl
610 /* Z = 8 */ {StableCluster, StableCl
611 },
612 { /* S = -2 */
613 /* A = 0
614 /* Z = 0 */ {StableCluster, StableCl
615 /* Z = 1 */ {StableCluster, StableCl
616 /* Z = 2 */ {StableCluster, StableCl
617 /* Z = 3 */ {StableCluster, StableCl
618 /* Z = 4 */ {StableCluster, StableCl
619 /* Z = 5 */ {StableCluster, StableCl
620 /* Z = 6 */ {StableCluster, StableCl
621 /* Z = 7 */ {StableCluster, StableCl
622 /* Z = 8 */ {StableCluster, StableCl
623 },
624 { /* S = -3 */
625 /* A = 0
626 /* Z = 0 */ {StableCluster, StableCl
627 /* Z = 1 */ {StableCluster, StableCl
628 /* Z = 2 */ {StableCluster, StableCl
629 /* Z = 3 */ {StableCluster, StableCl
630 /* Z = 4 */ {StableCluster, StableCl
631 /* Z = 5 */ {StableCluster, StableCl
632 /* Z = 6 */ {StableCluster, StableCl
633 /* Z = 7 */ {StableCluster, StableCl
634 /* Z = 8 */ {StableCluster, StableCl
635 }};
636
637 ParticleList decay(Cluster * const c) {
638 ParticleList decayProducts;
639 recursiveDecay(c, &decayProducts);
640
641 for(ParticleIter i = decayProducts.begin
642
643 // Correctly update the particle type
644 if(c->getA()==1) {
645 // assert(c->getZ()==1 || c->getZ()==0);
646 if(c->getZ()==1)
647 c->setType(Proton);
648 else if(c->getS()==-1)
649 c->setType(Lambda);
650 else
651 c->setType(Neutron);
652 c->setRealMass();
653 }
654
655 return decayProducts;
656 }
657
658 } // namespace ClusterDecay
659 } // namespace G4INCL
660
661