| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France << 27 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
28 // Joseph Cugnon, University of Liege, Belgium << 28 // Davide Mancusi, CEA
29 // Jean-Christophe David, CEA-Saclay, France << 29 // Alain Boudard, CEA
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H << 30 // Sylvie Leray, CEA
31 // Sylvie Leray, CEA-Saclay, France << 31 // Joseph Cugnon, University of Liege
32 // Davide Mancusi, CEA-Saclay, France << 32 //
>> 33 // INCL++ revision: v5.0_rc3
33 // 34 //
34 #define INCLXX_IN_GEANT4_MODE 1 35 #define INCLXX_IN_GEANT4_MODE 1
35 36
36 #include "globals.hh" 37 #include "globals.hh"
37 38
38 /** \file G4INCLClusterDecay.cc 39 /** \file G4INCLClusterDecay.cc
39 * \brief Static class for carrying out cluste 40 * \brief Static class for carrying out cluster decays
40 * 41 *
41 * \date 6th July 2011 << 42 * Created on: 6th July 2011
42 * \author Davide Mancusi << 43 * Author: Davide Mancusi
43 */ 44 */
44 45
45 #include "G4INCLClusterDecay.hh" 46 #include "G4INCLClusterDecay.hh"
46 #include "G4INCLParticleTable.hh" 47 #include "G4INCLParticleTable.hh"
47 #include "G4INCLKinematicsUtils.hh" 48 #include "G4INCLKinematicsUtils.hh"
48 #include "G4INCLRandom.hh" 49 #include "G4INCLRandom.hh"
49 #include "G4INCLPhaseSpaceGenerator.hh" << 50 //#include <cassert>
50 // #include <cassert> <<
51 #include <algorithm> <<
52 51
53 namespace G4INCL { 52 namespace G4INCL {
54 53
55 namespace ClusterDecay { << 54 ParticleList ClusterDecay::decay(Cluster * const c) {
56 << 55 ParticleList decayProducts;
57 namespace { << 56 recursiveDecay(c, &decayProducts);
58 << 57 return decayProducts;
59 /// \brief Carries out two-body decays << 58 }
60 void twoBodyDecay(Cluster * const c, Clu << 59
61 Particle *decayParticle = 0; << 60 void ClusterDecay::recursiveDecay(Cluster * const c, ParticleList *decayProducts) {
62 const ThreeVector mom(0.0, 0.0, 0.0); << 61 const G4int Z = c->getZ();
63 const ThreeVector pos = c->getPosition << 62 const G4int A = c->getA();
64 << 63
65 // Create the emitted particle << 64 ParticleTable::ClusterDecayType theDecayMode = ParticleTable::clusterDecayMode[Z][A];
66 switch(theDecayMode) { << 65
67 case ProtonDecay: << 66 switch(theDecayMode) {
68 decayParticle = new Particle(Proto << 67 default:
69 break; << 68 ERROR("Unrecognized cluster-decay mode: " << theDecayMode << std::endl
70 case NeutronDecay: << 69 << c->prG4int());
71 decayParticle = new Particle(Neutr << 70 case ParticleTable::StableCluster:
72 break; << 71 // For stable clusters, just return
73 case AlphaDecay: << 72 return;
74 decayParticle = new Cluster(2,4,0, << 73 break;
75 break; << 74 case ParticleTable::ProtonDecay:
76 case LambdaDecay: << 75 case ParticleTable::NeutronDecay:
77 decayParticle = new Particle(Lambd << 76 case ParticleTable::AlphaDecay:
78 break; << 77 // Two-body decays
79 default: << 78 twoBodyDecay(c, theDecayMode, decayProducts);
80 INCL_ERROR("Unrecognized cluster-d << 79 break;
81 << c->print()); << 80 case ParticleTable::TwoProtonDecay:
82 return; << 81 case ParticleTable::TwoNeutronDecay:
83 } << 82 // Three-body decays
84 decayParticle->makeParticipant(); << 83 threeBodyDecay(c, theDecayMode, decayProducts);
85 decayParticle->setNumberOfDecays(1); << 84 break;
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 } 85 }
576 86
577 /** \brief Table for cluster decays << 87 // Calls itself recursively in case the produced remnant is still unstable.
578 * << 88 // Sneaky, isn't it.
579 * Definition of "Stable": halflife > 1 ms << 89 recursiveDecay(c,decayProducts);
580 * << 90 }
581 * These table includes decay data for clu << 91
582 * not produce. It can't hurt. << 92 void ClusterDecay::twoBodyDecay(Cluster * const c, ParticleTable::ClusterDecayType theDecayMode, ParticleList *decayProducts) {
583 * << 93 Particle *decayParticle = 0;
584 * Unphysical nuclides (A<Z) are marked as << 94 const ThreeVector mom(0.0, 0.0, 0.0);
585 * produced by INCL. If you find them in t << 95 const ThreeVector pos = c->getPosition();
586 */ << 96
587 G4ThreadLocal ClusterDecayType clusterDeca << 97 // Create the emitted particle
588 {{/* S = 0 */ << 98 switch(theDecayMode) {
589 /* A = 0 << 99 case ParticleTable::ProtonDecay:
590 /* Z = 0 */ {StableCluster, StableCl << 100 decayParticle = new Particle(Proton, mom, pos);
591 /* Z = 1 */ {StableCluster, StableCl << 101 break;
592 /* Z = 2 */ {StableCluster, StableCl << 102 case ParticleTable::NeutronDecay:
593 /* Z = 3 */ {StableCluster, StableCl << 103 decayParticle = new Particle(Neutron, mom, pos);
594 /* Z = 4 */ {StableCluster, StableCl << 104 break;
595 /* Z = 5 */ {StableCluster, StableCl << 105 case ParticleTable::AlphaDecay:
596 /* Z = 6 */ {StableCluster, StableCl << 106 decayParticle = new Cluster(2,4);
597 /* Z = 7 */ {StableCluster, StableCl << 107 break;
598 /* Z = 8 */ {StableCluster, StableCl << 108 default:
599 }, << 109 ERROR("Unrecognized cluster-decay mode in two-body decay: " << theDecayMode << std::endl
600 { /* S = -1 */ << 110 << c->prG4int());
601 /* A = 0 << 111 return;
602 /* Z = 0 */ {StableCluster, StableCl << 112 }
603 /* Z = 1 */ {StableCluster, StableCl << 113 decayParticle->makeParticipant();
604 /* Z = 2 */ {StableCluster, StableCl << 114 decayParticle->setNumberOfDecays(1);
605 /* Z = 3 */ {StableCluster, StableCl << 115 decayParticle->setPosition(c->getPosition());
606 /* Z = 4 */ {StableCluster, StableCl << 116 decayParticle->setEmissionTime(c->getEmissionTime());
607 /* Z = 5 */ {StableCluster, StableCl << 117
608 /* Z = 6 */ {StableCluster, StableCl << 118 // Save some variables of the mother cluster
609 /* Z = 7 */ {StableCluster, StableCl << 119 const G4double motherMass = c->getMass();
610 /* Z = 8 */ {StableCluster, StableCl << 120 const ThreeVector velocity = -c->boostVector();
611 }, << 121
612 { /* S = -2 */ << 122 // Characteristics of the daughter particle
613 /* A = 0 << 123 const G4int daughterZ = c->getZ() - decayParticle->getZ();
614 /* Z = 0 */ {StableCluster, StableCl << 124 const G4int daughterA = c->getA() - decayParticle->getA();
615 /* Z = 1 */ {StableCluster, StableCl << 125 const G4double daughterMass = ParticleTable::getMass(daughterA,daughterZ);
616 /* Z = 2 */ {StableCluster, StableCl << 126
617 /* Z = 3 */ {StableCluster, StableCl << 127 // The mother cluster becomes the daughter
618 /* Z = 4 */ {StableCluster, StableCl << 128 c->setZ(daughterZ);
619 /* Z = 5 */ {StableCluster, StableCl << 129 c->setA(daughterA);
620 /* Z = 6 */ {StableCluster, StableCl << 130 c->setMass(daughterMass);
621 /* Z = 7 */ {StableCluster, StableCl << 131
622 /* Z = 8 */ {StableCluster, StableCl << 132 const G4double decayMass = decayParticle->getMass();
623 }, << 133 // assert(motherMass > daughterMass + decayMass); // Q-value should be >0
624 { /* S = -3 */ << 134
625 /* A = 0 << 135 // Decay kinematics in the mother rest frame
626 /* Z = 0 */ {StableCluster, StableCl << 136 const G4double pCM = KinematicsUtils::momentumInCM(motherMass, daughterMass, decayMass);
627 /* Z = 1 */ {StableCluster, StableCl << 137 const ThreeVector momentum = Random::normVector(pCM);
628 /* Z = 2 */ {StableCluster, StableCl << 138 c->setMomentum(momentum);
629 /* Z = 3 */ {StableCluster, StableCl << 139 c->adjustEnergyFromMomentum();
630 /* Z = 4 */ {StableCluster, StableCl << 140 decayParticle->setMomentum(-momentum);
631 /* Z = 5 */ {StableCluster, StableCl << 141 decayParticle->adjustEnergyFromMomentum();
632 /* Z = 6 */ {StableCluster, StableCl << 142
633 /* Z = 7 */ {StableCluster, StableCl << 143 // Boost to the lab frame
634 /* Z = 8 */ {StableCluster, StableCl << 144 decayParticle->boost(velocity);
635 }}; << 145 c->boost(velocity);
636 << 146
637 ParticleList decay(Cluster * const c) { << 147 // Add the decay particle to the list of decay products
638 ParticleList decayProducts; << 148 decayProducts->push_back(decayParticle);
639 recursiveDecay(c, &decayProducts); << 149 }
640 << 150
641 for(ParticleIter i = decayProducts.begin << 151 void ClusterDecay::threeBodyDecay(Cluster * const c, ParticleTable::ClusterDecayType theDecayMode, ParticleList *decayProducts) {
642 << 152 Particle *decayParticle1 = 0;
643 // Correctly update the particle type << 153 Particle *decayParticle2 = 0;
644 if(c->getA()==1) { << 154 const ThreeVector mom(0.0, 0.0, 0.0);
645 // assert(c->getZ()==1 || c->getZ()==0); << 155 const ThreeVector pos = c->getPosition();
646 if(c->getZ()==1) << 156
647 c->setType(Proton); << 157 // Create the emitted particles
648 else if(c->getS()==-1) << 158 switch(theDecayMode) {
649 c->setType(Lambda); << 159 case ParticleTable::TwoProtonDecay:
650 else << 160 decayParticle1 = new Particle(Proton, mom, pos);
651 c->setType(Neutron); << 161 decayParticle2 = new Particle(Proton, mom, pos);
652 c->setRealMass(); << 162 break;
653 } << 163 case ParticleTable::TwoNeutronDecay:
654 << 164 decayParticle1 = new Particle(Neutron, mom, pos);
655 return decayProducts; << 165 decayParticle2 = new Particle(Neutron, mom, pos);
>> 166 break;
>> 167 default:
>> 168 ERROR("Unrecognized cluster-decay mode in three-body decay: " << theDecayMode << std::endl
>> 169 << c->prG4int());
>> 170 return;
656 } 171 }
>> 172 decayParticle1->makeParticipant();
>> 173 decayParticle2->makeParticipant();
>> 174 decayParticle1->setNumberOfDecays(1);
>> 175 decayParticle2->setNumberOfDecays(1);
>> 176
>> 177 // Save some variables of the mother cluster
>> 178 const G4double motherMass = c->getMass();
>> 179 const ThreeVector velocity = -c->boostVector();
>> 180
>> 181 // Masses and charges of the daughter particle and of the decay products
>> 182 const G4int decayZ1 = decayParticle1->getZ();
>> 183 const G4int decayA1 = decayParticle1->getA();
>> 184 const G4int decayZ2 = decayParticle2->getZ();
>> 185 const G4int decayA2 = decayParticle2->getA();
>> 186 const G4int decayZ = decayZ1 + decayZ2;
>> 187 const G4int decayA = decayA1 + decayA2;
>> 188 const G4int daughterZ = c->getZ() - decayZ;
>> 189 const G4int daughterA = c->getA() - decayA;
>> 190 const G4double decayMass1 = decayParticle1->getMass();
>> 191 const G4double decayMass2 = decayParticle2->getMass();
>> 192 const G4double daughterMass = ParticleTable::getMass(daughterA,daughterZ);
>> 193
>> 194 // Q-values
>> 195 const G4double qValue = motherMass - daughterMass - decayMass1 - decayMass2;
>> 196 // assert(qValue > 0.); // Q-value should be >0
>> 197 const G4double qValueB = qValue * Random::shoot();
>> 198
>> 199 // The decay particles behave as if they had more mass until the second
>> 200 // decay
>> 201 const G4double decayMass = decayMass1 + decayMass2 + qValueB;
>> 202
>> 203 /* Stage A: mother --> daughter + (decay1+decay2) */
>> 204
>> 205 // The mother cluster becomes the daughter
>> 206 c->setZ(daughterZ);
>> 207 c->setA(daughterA);
>> 208 c->setMass(daughterMass);
>> 209
>> 210 // Decay kinematics in the mother rest frame
>> 211 const G4double pCMA = KinematicsUtils::momentumInCM(motherMass, daughterMass, decayMass);
>> 212 const ThreeVector momentumA = Random::normVector(pCMA);
>> 213 c->setMomentum(momentumA);
>> 214 c->adjustEnergyFromMomentum();
>> 215 const ThreeVector decayBoostVector = momentumA/std::sqrt(decayMass*decayMass + momentumA.mag2());
>> 216
>> 217 /* Stage B: (decay1+decay2) --> decay1 + decay2 */
>> 218
>> 219 // Decay kinematics in the (decay1+decay2) rest frame
>> 220 const G4double pCMB = KinematicsUtils::momentumInCM(decayMass, decayMass1, decayMass2);
>> 221 const ThreeVector momentumB = Random::normVector(pCMB);
>> 222 decayParticle1->setMomentum(momentumB);
>> 223 decayParticle2->setMomentum(-momentumB);
>> 224 decayParticle1->adjustEnergyFromMomentum();
>> 225 decayParticle2->adjustEnergyFromMomentum();
>> 226
>> 227 // Boost decay1 and decay2 to the Stage-A decay frame
>> 228 decayParticle1->boost(decayBoostVector);
>> 229 decayParticle2->boost(decayBoostVector);
>> 230
>> 231 // Boost all particles to the lab frame
>> 232 decayParticle1->boost(velocity);
>> 233 decayParticle2->boost(velocity);
>> 234 c->boost(velocity);
>> 235
>> 236 // Add the decay particles to the list of decay products
>> 237 decayProducts->push_back(decayParticle1);
>> 238 decayProducts->push_back(decayParticle2);
>> 239 }
657 240
658 } // namespace ClusterDecay << 241 }
659 } // namespace G4INCL <<
660 242
661 243