source: trunk/source/processes/hadronic/models/radioactive_decay/src/G4NuclearDecayChannel.cc@ 893

Last change on this file since 893 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 24.4 KB
Line 
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// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
27//
28// MODULES: G4NuclearDecayChannel.cc
29//
30// Version: 0.b.4
31// Date: 14/04/00
32// Author: F Lei & P R Truscott
33// Organisation: DERA UK
34// Customer: ESA/ESTEC, NOORDWIJK
35// Contract: 12115/96/JG/NL Work Order No. 3
36//
37// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
38//
39// CHANGE HISTORY
40// --------------
41//
42// 29 February 2000, P R Truscott, DERA UK
43// 0.b.3 release.
44//
45// 18 October 2002, F Lei
46// modified link metheds in DecayIt() to G4PhotoEvaporation() in order to
47// use the new Internal Coversion feature.
48// 13 April 2000, F Lei, DERA UK
49// Changes made are:
50// 1) Use PhotonEvaporation instead of DiscreteGammaDeexcitation
51// 2) verbose control
52//
53// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
54///////////////////////////////////////////////////////////////////////////////
55//
56#include "G4NuclearLevelManager.hh"
57#include "G4NuclearLevelStore.hh"
58#include "G4NuclearDecayChannel.hh"
59#include "G4DynamicParticle.hh"
60#include "G4DecayProducts.hh"
61#include "G4DecayTable.hh"
62#include "G4PhysicsLogVector.hh"
63#include "G4ParticleChangeForRadDecay.hh"
64#include "G4IonTable.hh"
65
66#include "G4BetaFermiFunction.hh"
67#include "G4PhotonEvaporation.hh"
68#include "G4AtomicDeexcitation.hh"
69
70
71const G4double G4NuclearDecayChannel:: pTolerance = 0.001;
72const G4double G4NuclearDecayChannel:: levelTolerance = 2.0*keV;
73//const G4bool G4NuclearDecayChannel:: FermiOn = true;
74
75///////////////////////////////////////////////////////////////////////////////
76//
77//
78// Constructor for one decay product (the nucleus).
79//
80G4NuclearDecayChannel::G4NuclearDecayChannel
81 (const G4RadioactiveDecayMode &theMode,
82 G4int Verbose,
83 const G4ParticleDefinition *theParentNucleus,
84 G4double theBR,
85 G4double theQtransition,
86 G4int A,
87 G4int Z,
88 G4double theDaughterExcitation) :
89 G4GeneralPhaseSpaceDecay(Verbose), decayMode(theMode)
90{
91#ifdef G4VERBOSE
92 if (GetVerboseLevel()>1)
93 {G4cout <<"G4NuclearDecayChannel constructor for " <<G4int(theMode) <<G4endl;}
94#endif
95 SetParent(theParentNucleus);
96 FillParent();
97 parent_mass = theParentNucleus->GetPDGMass();
98 SetBR (theBR);
99 SetNumberOfDaughters (1);
100 FillDaughterNucleus (0, A, Z, theDaughterExcitation);
101 Qtransition = theQtransition;
102}
103///////////////////////////////////////////////////////////////////////////////
104//
105//
106// Constructor for a daughter nucleus and one other particle.
107//
108G4NuclearDecayChannel::G4NuclearDecayChannel
109 (const G4RadioactiveDecayMode &theMode,
110 G4int Verbose,
111 const G4ParticleDefinition *theParentNucleus,
112 G4double theBR,
113 G4double theQtransition,
114 G4int A,
115 G4int Z,
116 G4double theDaughterExcitation,
117 const G4String theDaughterName1) :
118 G4GeneralPhaseSpaceDecay(Verbose), decayMode(theMode)
119{
120#ifdef G4VERBOSE
121 if (GetVerboseLevel()>1)
122 {G4cout <<"G4NuclearDecayChannel constructor for " <<G4int(theMode) <<G4endl;}
123#endif
124 SetParent (theParentNucleus);
125 FillParent();
126 parent_mass = theParentNucleus->GetPDGMass();
127 SetBR (theBR);
128 SetNumberOfDaughters (2);
129 SetDaughter(0, theDaughterName1);
130 FillDaughterNucleus (1, A, Z, theDaughterExcitation);
131 Qtransition = theQtransition;
132}
133///////////////////////////////////////////////////////////////////////////////
134//
135//
136// Constructor for a daughter nucleus and two other particles.
137//
138G4NuclearDecayChannel::G4NuclearDecayChannel
139 (const G4RadioactiveDecayMode &theMode,
140 G4int Verbose,
141 const G4ParticleDefinition *theParentNucleus,
142 G4double theBR,
143 G4double theFFN,
144 G4bool betaS,
145 CLHEP::RandGeneral* randBeta,
146 G4double theQtransition,
147 G4int A,
148 G4int Z,
149 G4double theDaughterExcitation,
150 const G4String theDaughterName1,
151 const G4String theDaughterName2) :
152 G4GeneralPhaseSpaceDecay(Verbose), decayMode(theMode)
153 //,BetaSimple(betaS),
154 // RandomEnergy(randBeta), Qtransition(theQtransition),FermiFN(theFFN)
155{
156#ifdef G4VERBOSE
157 if (GetVerboseLevel()>1)
158 {G4cout <<"G4NuclearDecayChannel constructor for " <<G4int(theMode) <<G4endl;}
159#endif
160 SetParent (theParentNucleus);
161 FillParent();
162 parent_mass = theParentNucleus->GetPDGMass();
163 SetBR (theBR);
164 SetNumberOfDaughters (3);
165 SetDaughter(0, theDaughterName1);
166 SetDaughter(2, theDaughterName2);
167 FillDaughterNucleus(1, A, Z, theDaughterExcitation);
168 BetaSimple = betaS;
169 RandomEnergy = randBeta;
170 Qtransition = theQtransition;
171 FermiFN = theFFN;
172}
173
174////////////////////////////////////////////////////////////////////////////////
175//
176//
177//
178//
179#include "G4HadTmpUtil.hh"
180
181void G4NuclearDecayChannel::FillDaughterNucleus (G4int index, G4int A, G4int Z,
182 G4double theDaughterExcitation)
183{
184 //
185 //
186 // Determine if the proposed daughter nucleus has a sensible A, Z and excitation
187 // energy.
188 //
189 if (A<1 || Z<0 || theDaughterExcitation <0.0)
190 {
191 G4cerr <<"Error in G4NuclearDecayChannel::FillDaughterNucleus";
192 G4cerr <<"Inappropriate values of daughter A, Z or excitation" <<G4endl;
193 G4cerr <<"A = " <<A <<" and Z = " <<Z;
194 G4cerr <<" Ex = " <<theDaughterExcitation*MeV <<"MeV" <<G4endl;
195 G4Exception(__FILE__, G4inttostring(__LINE__), FatalException, "G4NuclearDecayChannel::FillDaughterNucleus");
196 }
197 //
198 //
199 // Save A and Z to local variables. Find the GROUND STATE of the daughter
200 // nucleus and save this, as an ion, in the array of daughters.
201 //
202 daughterA = A;
203 daughterZ = Z;
204 if (Z == 1 && A == 1) {
205 daughterNucleus = G4Proton::Definition();
206 } else if (Z == 0 && A == 1) {
207 daughterNucleus = G4Neutron::Definition();
208 } else {
209 G4IonTable *theIonTable =
210 (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable());
211 daughterNucleus = theIonTable->GetIon(daughterZ, daughterA, theDaughterExcitation*MeV);
212 }
213 daughterExcitation = theDaughterExcitation;
214 SetDaughter(index, daughterNucleus);
215}
216
217///////////////////////////////////////////////////////////////////////////////
218//
219//
220//
221//
222G4DecayProducts *G4NuclearDecayChannel::DecayIt (G4double theParentMass)
223{
224 //
225 //
226 // Load-up the details of the parent and daughter particles if they have not
227 // been defined properly.
228 //
229 if (parent == 0) FillParent();
230 if (daughters == 0) FillDaughters();
231 //
232 //
233 // We want to ensure that the difference between the total
234 // parent and daughter masses equals the energy liberated by the transition.
235 //
236 theParentMass = 0.0;
237 for( G4int index=0; index < numberOfDaughters; index++)
238 {theParentMass += daughters[index]->GetPDGMass();}
239 theParentMass += Qtransition ;
240 // bug fix for beta+ decay (flei 25/09/01)
241 if (decayMode == 2) theParentMass -= 2*0.511 * MeV;
242 //
243 if (GetVerboseLevel()>1) {
244 G4cout << "G4NuclearDecayChannel::DecayIt ";
245 G4cout << "the decay mass = " << theParentMass << G4endl;
246 }
247
248 SetParentMass (theParentMass);
249
250 //
251 //
252 // Define a product vector.
253 //
254 G4DecayProducts *products = 0;
255 //
256 //
257 // Depending upon the number of daughters, select the appropriate decay
258 // kinematics scheme.
259 //
260 switch (numberOfDaughters)
261 {
262 case 0:
263 if (GetVerboseLevel()>0)
264 {
265 G4cout << "G4NuclearDecayChannel::DecayIt ";
266 G4cout << " daughters not defined " <<G4endl;
267 }
268 break;
269 case 1:
270 products = OneBodyDecayIt();
271 break;
272 case 2:
273 products = TwoBodyDecayIt();
274 break;
275 case 3:
276 products = BetaDecayIt();
277 break;
278 default:
279 G4cerr <<"Error in G4NuclearDecayChannel::DecayIt" <<G4endl;
280 G4cerr <<"Number of daughters in decay = " <<numberOfDaughters <<G4endl;
281 G4Exception(__FILE__, G4inttostring(__LINE__), FatalException, "G4NuclearDecayChannel::DecayIt");
282 }
283 if ((products == 0) && (GetVerboseLevel()>0)) {
284 G4cerr << "G4NuclearDecayChannel::DecayIt ";
285 G4cerr << *parent_name << " can not decay " << G4endl;
286 DumpInfo();
287 }
288
289 //
290 // now we have to take care of the EC product which have go through the ARM
291 if (decayMode == 3 || decayMode == 4 || decayMode == 5) {
292 G4int eShell = 0;
293 switch (decayMode)
294 {
295 case KshellEC:
296 //
297 {
298 eShell = 1;
299 }
300 break;
301 case LshellEC:
302 //
303 {
304 eShell = G4int(G4UniformRand()*3)+1;
305 }
306 break;
307 case MshellEC:
308 //
309 {
310 eShell = G4int(G4UniformRand()*5)+4;
311 }
312 break;
313 case ERROR:
314 default:
315 G4cout << " There is an error in decay mode selection! exit RDM now" << G4endl;
316 exit(0);
317 }
318 G4int aZ = daughterZ;
319 if (aZ > 5 && aZ < 101) { // only applies to 5< Z <101
320 G4AtomicDeexcitation* atomDeex = new G4AtomicDeexcitation();
321 //no Auger electron generation
322 // atomDeex->ActivateAugerElectronProduction(0);
323 std::vector<G4DynamicParticle*>* armProducts = atomDeex->GenerateParticles(aZ,eShell);
324
325 // pop up the daughter before insertion
326 dynamicDaughter = products->PopProducts();
327 for (size_t i = 0; i < armProducts->size(); i++)
328 products->PushProducts ((*armProducts)[i]);
329 delete armProducts;
330 delete atomDeex;
331 products->PushProducts (dynamicDaughter);
332 }
333 }
334
335 //
336 // If the decay is to an excited state of the daughter nuclide, we need
337 // to apply the photo-evaporation process.
338 //
339 if (daughterExcitation > 0.0)
340 {
341 //
342 //
343 // Pop the daughter nucleus off the product vector - we need to retain
344 // the momentum of this particle.
345 //
346 dynamicDaughter = products->PopProducts();
347 G4LorentzVector daughterMomentum = dynamicDaughter->Get4Momentum();
348 G4ThreeVector const daughterMomentum1(static_cast<const G4LorentzVector> (daughterMomentum));
349 //
350 //
351 // Now define a G4Fragment with the correct A, Z and excitation, and declare and
352 // initialise a G4DiscreteGammaDeexcitation object.
353 //
354
355 // daughterMomentum.setT(daughterMomentum.t()+G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass( daughterZ, daughterA )+daughterExcitation);
356
357 // daughterMomentum.setT(daughterMomentum.t()+daughterExcitation);
358 G4Fragment nucleus(daughterA, daughterZ, daughterMomentum);
359 //G4LorentzVector p4(0.,0.,0.,G4NucleiProperties::GetNuclearMass(daughterA,daughterZ)
360 // +daughterExcitation);
361 //G4Fragment nucleus(daughterA, daughterZ, p4);
362 // nucleus.SetExcitationEnergy(daughterExcitation);
363
364 // G4VGammaDeexcitation* deexcitation = new G4DiscreteGammaDeexcitation;
365 G4PhotonEvaporation* deexcitation = new G4PhotonEvaporation;
366 deexcitation->SetVerboseLevel(GetVerboseLevel());
367 // deexcitation->Initialize(nucleus);
368 deexcitation->SetICM(true);
369 if (decayMode == 0) {
370 deexcitation->RDMForced(true);
371 } else {
372 deexcitation->RDMForced(false);
373 }
374 // ARM in G4 is applied but no auger electrons!
375 deexcitation->SetARM(true);
376 // not applied
377 //deexcitation->SetARM(false);
378 //
379 deexcitation->SetMaxHalfLife(1e-6*second);
380 //
381 // Get the gammas by deexciting the nucleus.
382 //
383 G4FragmentVector* gammas = deexcitation->BreakItUp(nucleus);
384 // in the case of BreakItUp(nucleus), the returned G4FragmentVector contains the residual nuclide
385 // as its last entry.
386 G4int nGammas=gammas->size()-1;
387 //
388 //
389 // Go through each gamma/e- and add it to the decay product. The angular distribution
390 // of the gammas is isotropic, and the residual nucleus is assumed not to suffer
391 // any recoil as a result of this de-excitation.
392 //
393 for (G4int ig=0; ig<nGammas; ig++)
394 {
395 // G4double costheta = 2.0*G4UniformRand() - 1.0;
396 // G4double sintheta = std::sqrt((1.0 - costheta) * (1.0+costheta));
397 // G4double phi = twopi * G4UniformRand();
398 // G4ParticleMomentum gDirection
399 // (sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
400 //G4double gEnergy = gammas->operator[](ig)->GetMomentum().e()
401 // - gammas->operator[](ig)->GetParticleDefinition()->GetPDGMass() ;
402 G4DynamicParticle *theGammaRay = new
403 G4DynamicParticle (gammas->operator[](ig)->GetParticleDefinition(),
404 gammas->operator[](ig)->GetMomentum());
405 theGammaRay -> SetProperTime(gammas->operator[](ig)->GetCreationTime());
406 products->PushProducts (theGammaRay);
407 }
408 //
409 // now the nucleus
410 G4double finalDaughterExcitation = gammas->operator[](nGammas)->GetExcitationEnergy();
411 // f.lei (03/01/03) this is needed to fix the crach in test18
412 if (finalDaughterExcitation <= 1.0*keV) finalDaughterExcitation = 0 ;
413 G4IonTable *theIonTable = (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable());
414 // f.lei (07/03/05) added the delete to fix bug#711
415 if (dynamicDaughter) delete dynamicDaughter;
416 dynamicDaughter = new G4DynamicParticle
417 (theIonTable->GetIon(daughterZ,daughterA,finalDaughterExcitation),
418 daughterMomentum1);
419 products->PushProducts (dynamicDaughter);
420 //
421 // Delete/reset variables associated with the gammas.
422 //
423 // if (nGammas != 0) gammas->clearAndDestroy();
424 while (!gammas->empty()) {
425 delete *(gammas->end()-1);
426 gammas->pop_back();
427 }
428 // gammas->clearAndDestroy();
429 delete gammas;
430 delete deexcitation;
431 }
432 return products;
433}
434////////////////////////////////////////////////////////////////////////////////
435//
436
437G4DecayProducts *G4NuclearDecayChannel::BetaDecayIt()
438
439{
440 if (GetVerboseLevel()>1) G4cout << "G4Decay::BetaDecayIt()"<<G4endl;
441
442 //daughters'mass
443 G4double daughtermass[3];
444 G4double sumofdaughtermass = 0.0;
445 G4double pmass = GetParentMass();
446 for (G4int index=0; index<3; index++)
447 {
448 daughtermass[index] = daughters[index]->GetPDGMass();
449 sumofdaughtermass += daughtermass[index];
450 }
451
452 //create parent G4DynamicParticle at rest
453 G4ParticleMomentum dummy;
454 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
455
456 //create G4Decayproducts
457 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
458 delete parentparticle;
459
460 G4double Q = pmass - sumofdaughtermass;
461
462 // 09/11/2004 flei
463 // All Beta decays are now treated with the improved 3 body decay algorithm. No more slow/fast modes
464 /*
465 if (BetaSimple == true) {
466
467 // Use the histogramed distribution to generate the beta energy
468 G4double daughtermomentum[2];
469 G4double daughterenergy[2];
470 daughterenergy[0] = RandomEnergy->shoot() * Q;
471 daughtermomentum[0] = std::sqrt(daughterenergy[0]*daughterenergy[0] +
472 2.0*daughterenergy[0] * daughtermass[0]);
473 // the recoil neuleus is asummed to have a maximum energy of Q/daughterA/1000.
474 daughterenergy[1] = G4UniformRand() * Q/(1000.*daughterA);
475 G4double recoilmomentumsquared = daughterenergy[1]*daughterenergy[1] +
476 2.0*daughterenergy[1] * daughtermass[1];
477 if (recoilmomentumsquared < 0.0) recoilmomentumsquared = 0.0;
478 daughtermomentum[1] = std::sqrt(recoilmomentumsquared);
479 //
480 //create daughter G4DynamicParticle
481 G4double costheta, sintheta, phi, sinphi, cosphi;
482 // G4double costhetan, sinthetan, phin, sinphin, cosphin;
483 costheta = 2.*G4UniformRand()-1.0;
484 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
485 phi = twopi*G4UniformRand()*rad;
486 sinphi = std::sin(phi);
487 cosphi = std::cos(phi);
488 G4ParticleMomentum direction0(sintheta*cosphi,sintheta*sinphi,costheta);
489 G4DynamicParticle * daughterparticle
490 = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
491 products->PushProducts(daughterparticle);
492 // The two products are independent in directions
493 costheta = 2.*G4UniformRand()-1.0;
494 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
495 phi = twopi*G4UniformRand()*rad;
496 sinphi = std::sin(phi);
497 cosphi = std::cos(phi);
498 G4ParticleMomentum direction1(sintheta*cosphi,sintheta*sinphi,costheta);
499 daughterparticle
500 = new G4DynamicParticle( daughters[1], direction1*daughtermomentum[1]);
501 products->PushProducts(daughterparticle);
502
503 // the neutrino is igored in this case
504
505 } else {
506 */
507 /* original slow method
508 //calculate daughter momentum
509 // Generate two
510 G4double rd1, rd2;
511 G4double daughtermomentum[3];
512 G4double daughterenergy[3];
513 G4double momentummax=0.0, momentumsum = 0.0;
514 G4double fermif;
515 G4BetaFermiFunction* aBetaFermiFunction;
516 if (decayMode == 1) {
517 // beta-decay
518 aBetaFermiFunction = new G4BetaFermiFunction (daughterA, daughterZ);
519 } else {
520 // beta+decay
521 aBetaFermiFunction = new G4BetaFermiFunction (daughterA, -daughterZ);
522 }
523 if (GetVerboseLevel()>1) {
524 G4cout<< " Q = " <<Q<<G4endl;
525 G4cout<< " daughterA = " <<daughterA<<G4endl;
526 G4cout<< " daughterZ = " <<daughterZ<<G4endl;
527 G4cout<< " decayMode = " <<static_cast<G4int>(decayMode) << G4endl;
528 G4cout<< " FermiFN = " <<FermiFN<<G4endl;
529 }
530 do
531 {
532 rd1 = G4UniformRand();
533 rd2 = G4UniformRand();
534
535 momentummax = 0.0;
536 momentumsum = 0.0;
537
538 // daughter 0
539
540 // energy = rd2*(pmass - sumofdaughtermass);
541 daughtermomentum[0] = std::sqrt(rd2) * std::sqrt((Q + 2.0*daughtermass[0])*Q);
542 daughterenergy[0] = std::sqrt(daughtermomentum[0]*daughtermomentum[0] +
543 daughtermass[0] * daughtermass[0]) - daughtermass[0];
544 if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0];
545 momentumsum += daughtermomentum[0];
546
547 // daughter 2
548 // energy = (1.-rd1)*(pmass - sumofdaughtermass);
549 daughtermomentum[2] = std::sqrt(rd1)*std::sqrt((Q + 2.0*daughtermass[2])*Q);
550 daughterenergy[2] = std::sqrt(daughtermomentum[2]*daughtermomentum[2] +
551 daughtermass[2] * daughtermass[2]) - daughtermass[2];
552 if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2];
553 momentumsum += daughtermomentum[2];
554
555 // daughter 1
556
557 daughterenergy[1] = Q - daughterenergy[0] - daughterenergy[2];
558 if (daughterenergy[1] > 0.0) {
559 daughtermomentum[1] = std::sqrt(daughterenergy[1]*daughterenergy[1] +
560 2.0*daughterenergy[1] * daughtermass[1]);
561 if ( daughtermomentum[1] >momentummax ) momentummax =
562 daughtermomentum[1];
563 momentumsum += daughtermomentum[1];
564 } else {
565 momentummax = momentumsum = Q;
566 }
567 // beta particles is sampled with no coulomb effects applied above. Now
568 // apply the Fermi function using rejection method.
569 daughterenergy[0] = daughterenergy[0]*MeV/0.511;
570 fermif = aBetaFermiFunction->GetFF(daughterenergy[0])/FermiFN;
571 // fermif: normalised Fermi factor
572 if (G4UniformRand() > fermif) momentummax = momentumsum = Q;
573 // rejection method
574 } while (momentummax > momentumsum - momentummax );
575 delete aBetaFermiFunction;
576
577 // output message
578 if (GetVerboseLevel()>1) {
579 G4cout <<" daughter 0:" <<daughtermomentum[0]/GeV <<"[GeV/c]" <<G4endl;
580 G4cout <<" daughter 1:" <<daughtermomentum[1]/GeV <<"[GeV/c]" <<G4endl;
581 G4cout <<" daughter 2:" <<daughtermomentum[2]/GeV <<"[GeV/c]" <<G4endl;
582 G4cout <<" momentum sum:" <<momentumsum/GeV <<"[GeV/c]" <<G4endl;
583 }
584 */
585 // faster method as suggested by Dirk Kruecker of FZ-Julich
586 G4double daughtermomentum[3];
587 G4double daughterenergy[3];
588 // Use the histogramed distribution to generate the beta energy
589 daughterenergy[0] = RandomEnergy->shoot() * Q;
590 daughtermomentum[0] = std::sqrt(daughterenergy[0]*daughterenergy[0] +
591 2.0*daughterenergy[0] * daughtermass[0]);
592
593 //neutrino energy distribution is flat within the kinematical limits
594 G4double rd = 2*G4UniformRand()-1;
595 // limits
596 G4double Mme=pmass-daughtermass[0];
597 G4double K=0.5-daughtermass[1]*daughtermass[1]/(2*Mme*Mme-4*pmass*daughterenergy[0]);
598
599 daughterenergy[2]=K*(Mme-daughterenergy[0]+rd*daughtermomentum[0]);
600 daughtermomentum[2] = daughterenergy[2] ;
601
602 // the recoil neuleus
603 daughterenergy[1] = Q-daughterenergy[0]-daughterenergy[2];
604 G4double recoilmomentumsquared = daughterenergy[1]*daughterenergy[1] +
605 2.0*daughterenergy[1] * daughtermass[1];
606 if (recoilmomentumsquared < 0.0) recoilmomentumsquared = 0.0;
607 daughtermomentum[1] = std::sqrt(recoilmomentumsquared);
608
609 // output message
610 if (GetVerboseLevel()>1) {
611 G4cout <<" daughter 0:" <<daughtermomentum[0]/GeV <<"[GeV/c]" <<G4endl;
612 G4cout <<" daughter 1:" <<daughtermomentum[1]/GeV <<"[GeV/c]" <<G4endl;
613 G4cout <<" daughter 2:" <<daughtermomentum[2]/GeV <<"[GeV/c]" <<G4endl;
614 }
615 //create daughter G4DynamicParticle
616 G4double costheta, sintheta, phi, sinphi, cosphi;
617 G4double costhetan, sinthetan, phin, sinphin, cosphin;
618 costheta = 2.*G4UniformRand()-1.0;
619 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
620 phi = twopi*G4UniformRand()*rad;
621 sinphi = std::sin(phi);
622 cosphi = std::cos(phi);
623 G4ParticleMomentum direction0(sintheta*cosphi,sintheta*sinphi,costheta);
624 G4DynamicParticle * daughterparticle
625 = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
626 products->PushProducts(daughterparticle);
627
628 costhetan = (daughtermomentum[1]*daughtermomentum[1]-
629 daughtermomentum[2]*daughtermomentum[2]-
630 daughtermomentum[0]*daughtermomentum[0])/
631 (2.0*daughtermomentum[2]*daughtermomentum[0]);
632 // added the following test to avoid rounding erros. A problem
633 // reported bye Ben Morgan of Uni.Warwick
634 if (costhetan > 1.) costhetan = 1.;
635 if (costhetan < -1.) costhetan = -1.;
636 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
637 phin = twopi*G4UniformRand()*rad;
638 sinphin = std::sin(phin);
639 cosphin = std::cos(phin);
640 G4ParticleMomentum direction2;
641 direction2.setX( sinthetan*cosphin*costheta*cosphi -
642 sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
643 direction2.setY( sinthetan*cosphin*costheta*sinphi +
644 sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
645 direction2.setZ( -sinthetan*cosphin*sintheta +
646 costhetan*costheta);
647 daughterparticle = new G4DynamicParticle
648 ( daughters[2], direction2*(daughtermomentum[2]/direction2.mag()));
649 products->PushProducts(daughterparticle);
650
651 daughterparticle =
652 new G4DynamicParticle (daughters[1],
653 (direction0*daughtermomentum[0] +
654 direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0));
655 products->PushProducts(daughterparticle);
656 // }
657 // delete daughterparticle;
658
659 if (GetVerboseLevel()>1) {
660 G4cout << "G4NuclearDecayChannel::BetaDecayIt ";
661 G4cout << " create decay products in rest frame " <<G4endl;
662 products->DumpInfo();
663 }
664 return products;
665}
Note: See TracBrowser for help on using the repository browser.