Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4PhaseSpaceDecayChannel.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//
27// $Id$
28//
29//
30// ------------------------------------------------------------
31// GEANT 4 class header file
32//
33// History: first implementation, based on object model of
34// 27 July 1996 H.Kurashige
35// 20 Oct 1996 H.Kurashige
36// 30 May 1997 H.Kurashige
37// ------------------------------------------------------------
38
41#include "G4SystemOfUnits.hh"
42#include "G4DecayProducts.hh"
43#include "G4VDecayChannel.hh"
45#include "Randomize.hh"
46#include "G4LorentzVector.hh"
47#include "G4LorentzRotation.hh"
48
49
51 :G4VDecayChannel("Phase Space", Verbose)
52{
53
54}
55
57 const G4String& theParentName,
58 G4double theBR,
59 G4int theNumberOfDaughters,
60 const G4String& theDaughterName1,
61 const G4String& theDaughterName2,
62 const G4String& theDaughterName3,
63 const G4String& theDaughterName4 )
64 :G4VDecayChannel("Phase Space",
65 theParentName,theBR,
66 theNumberOfDaughters,
67 theDaughterName1,
68 theDaughterName2,
69 theDaughterName3,
70 theDaughterName4),
71 current_parent_mass(0.0)
72{
73
74}
75
77{
78}
79
81{
82#ifdef G4VERBOSE
83 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
84#endif
85
86 G4DecayProducts * products = 0;
87
88 if (parent == 0) FillParent();
89 if (daughters == 0) FillDaughters();
90
91 if (parentMass >0.0) current_parent_mass = parentMass;
92 else current_parent_mass = parent_mass;
93
94 switch (numberOfDaughters){
95 case 0:
96#ifdef G4VERBOSE
97 if (GetVerboseLevel()>0) {
98 G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
99 G4cout << " daughters not defined " <<G4endl;
100 }
101#endif
102 break;
103 case 1:
104 products = OneBodyDecayIt();
105 break;
106 case 2:
107 products = TwoBodyDecayIt();
108 break;
109 case 3:
110 products = ThreeBodyDecayIt();
111 break;
112 default:
113 products = ManyBodyDecayIt();
114 break;
115 }
116#ifdef G4VERBOSE
117 if ((products == 0) && (GetVerboseLevel()>0)) {
118 G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
119 G4cout << *parent_name << " can not decay " << G4endl;
120 DumpInfo();
121 }
122#endif
123 return products;
124}
125
126G4DecayProducts *G4PhaseSpaceDecayChannel::OneBodyDecayIt()
127{
128#ifdef G4VERBOSE
129 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt()"<<G4endl;
130#endif
131
132 //create parent G4DynamicParticle at rest
133 G4ThreeVector dummy;
134 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
135 //create G4Decayproducts
136 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
137 delete parentparticle;
138
139 //create daughter G4DynamicParticle at rest
140 G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0], dummy, 0.0);
141 products->PushProducts(daughterparticle);
142
143#ifdef G4VERBOSE
144 if (GetVerboseLevel()>1) {
145 G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt ";
146 G4cout << " create decay products in rest frame " <<G4endl;
147 products->DumpInfo();
148 }
149#endif
150 return products;
151}
152
153G4DecayProducts *G4PhaseSpaceDecayChannel::TwoBodyDecayIt()
154{
155#ifdef G4VERBOSE
156 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt()"<<G4endl;
157#endif
158 // parent mass
159 G4double parentmass = current_parent_mass;
160
161 //daughters'mass
162 G4double daughtermass[2];
163 G4double daughtermomentum;
164 daughtermass[0] = daughters_mass[0];
165 daughtermass[1] = daughters_mass[1];
166
167 //create parent G4DynamicParticle at rest
168 G4ThreeVector dummy;
169 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
170 //create G4Decayproducts
171 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
172 delete parentparticle;
173
174 //calculate daughter momentum
175 daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
176 if (daughtermomentum <0.0) {
177#ifdef G4VERBOSE
178 if (GetVerboseLevel()>0) {
179 G4cerr << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt "
180 << "sum of daughter mass is larger than parent mass" << G4endl;
181 G4cerr << "parent :" << parent->GetParticleName() << " " << current_parent_mass/GeV << G4endl;
182 G4cerr << "daughter 1 :" << daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
183 G4cerr << "daughter 2:" << daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
184 }
185#endif
186 G4Exception("G4PhaseSpaceDecayChannel::TwoBodyDecayIt",
187 "PART112", JustWarning,
188 "Can not create decay products: sum of daughter mass is larger than parent mass");
189 return products;
190 }
191
192 G4double costheta = 2.*G4UniformRand()-1.0;
193 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
194 G4double phi = twopi*G4UniformRand()*rad;
195 G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
196
197 //create daughter G4DynamicParticle
198 G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0], direction*daughtermomentum);
199 products->PushProducts(daughterparticle);
200 daughterparticle = new G4DynamicParticle( daughters[1], direction*(-1.0*daughtermomentum));
201 products->PushProducts(daughterparticle);
202
203#ifdef G4VERBOSE
204 if (GetVerboseLevel()>1) {
205 G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt ";
206 G4cout << " create decay products in rest frame " <<G4endl;
207 products->DumpInfo();
208 }
209#endif
210 return products;
211}
212
213G4DecayProducts *G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()
214// algorism of this code is originally written in GDECA3 of GEANT3
215{
216#ifdef G4VERBOSE
217 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()"<<G4endl;
218#endif
219 // parent mass
220 G4double parentmass = current_parent_mass;
221 //daughters'mass
222 G4double daughtermass[3];
223 G4double sumofdaughtermass = 0.0;
224 for (G4int index=0; index<3; index++){
225 daughtermass[index] = daughters_mass[index];
226 sumofdaughtermass += daughtermass[index];
227 }
228
229 //create parent G4DynamicParticle at rest
230 G4ThreeVector dummy;
231 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
232
233
234 //create G4Decayproducts
235 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
236 delete parentparticle;
237
238 if (sumofdaughtermass >parentmass) {
239#ifdef G4VERBOSE
240 if (GetVerboseLevel()>0) {
241 G4cerr << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt "
242 << "sum of daughter mass is larger than parent mass" << G4endl;
243 G4cerr << "parent :" << parent->GetParticleName() << " " << current_parent_mass/GeV << G4endl;
244 G4cerr << "daughter 1 :" << daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
245 G4cerr << "daughter 2:" << daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
246 G4cerr << "daughter 3:" << daughters[2]->GetParticleName() << " " << daughtermass[2]/GeV << G4endl;
247 }
248#endif
249 G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
250 "PART112", JustWarning,
251 "Can not create decay products: sum of daughter mass is larger than parent mass");
252 return products;
253 }
254
255
256
257 //calculate daughter momentum
258 // Generate two
259 G4double rd1, rd2, rd;
260 G4double daughtermomentum[3];
261 G4double momentummax=0.0, momentumsum = 0.0;
262 G4double energy;
263
264 do {
265 rd1 = G4UniformRand();
266 rd2 = G4UniformRand();
267 if (rd2 > rd1) {
268 rd = rd1;
269 rd1 = rd2;
270 rd2 = rd;
271 }
272 momentummax = 0.0;
273 momentumsum = 0.0;
274 // daughter 0
275 energy = rd2*(parentmass - sumofdaughtermass);
276 daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
277 if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0];
278 momentumsum += daughtermomentum[0];
279 // daughter 1
280 energy = (1.-rd1)*(parentmass - sumofdaughtermass);
281 daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
282 if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1];
283 momentumsum += daughtermomentum[1];
284 // daughter 2
285 energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
286 daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
287 if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2];
288 momentumsum += daughtermomentum[2];
289 } while (momentummax > momentumsum - momentummax );
290 // output message
291#ifdef G4VERBOSE
292 if (GetVerboseLevel()>1) {
293 G4cout << " daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl;
294 G4cout << " daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl;
295 G4cout << " daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl;
296 G4cout << " momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl;
297 }
298#endif
299
300 //create daughter G4DynamicParticle
301 G4double costheta, sintheta, phi, sinphi, cosphi;
302 G4double costhetan, sinthetan, phin, sinphin, cosphin;
303 costheta = 2.*G4UniformRand()-1.0;
304 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
305 phi = twopi*G4UniformRand()*rad;
306 sinphi = std::sin(phi);
307 cosphi = std::cos(phi);
308 G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
309 G4DynamicParticle * daughterparticle
310 = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
311 products->PushProducts(daughterparticle);
312
313 costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
314 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
315 phin = twopi*G4UniformRand()*rad;
316 sinphin = std::sin(phin);
317 cosphin = std::cos(phin);
318 G4ThreeVector direction2;
319 direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
320 direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
321 direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
322 daughterparticle = new G4DynamicParticle( daughters[2], direction2*(daughtermomentum[2]/direction2.mag()));
323 products->PushProducts(daughterparticle);
324
325 daughterparticle =
327 daughters[1],
328 (direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0)
329 );
330 products->PushProducts(daughterparticle);
331
332#ifdef G4VERBOSE
333 if (GetVerboseLevel()>1) {
334 G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt ";
335 G4cout << " create decay products in rest frame " <<G4endl;
336 products->DumpInfo();
337 }
338#endif
339 return products;
340}
341
342G4DecayProducts *G4PhaseSpaceDecayChannel::ManyBodyDecayIt()
343// algorithm of this code is originally written in FORTRAN by M.Asai
344//******************************************************************
345// NBODY
346// N-body phase space Monte-Carlo generator
347// Makoto Asai
348// Hiroshima Institute of Technology
350// Revised release : 19/Apr/1995
351//
352{
353 G4int index, index2;
354
355 //return value
356 G4DecayProducts *products;
357
358#ifdef G4VERBOSE
359 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()"<<G4endl;
360#endif
361 // parent mass
362 G4double parentmass = current_parent_mass;
363 //daughters'mass
364 G4double *daughtermass = new G4double[numberOfDaughters];
365 G4double sumofdaughtermass = 0.0;
366 for (index=0; index<numberOfDaughters; index++){
367 daughtermass[index] = daughters_mass[index];
368 sumofdaughtermass += daughtermass[index];
369 }
370
371 //Calculate daughter momentum
372 G4double *daughtermomentum = new G4double[numberOfDaughters];
373 G4ThreeVector direction;
374 G4DynamicParticle **daughterparticle;
376 G4double tmas;
377 G4double weight = 1.0;
378 G4int numberOfTry = 0;
379
380 do {
381 //Generate rundom number in descending order
382 G4double temp;
384 rd[0] = 1.0;
385 for(index =1; index < numberOfDaughters -1; index++) rd[index] = G4UniformRand();
386 rd[ numberOfDaughters -1] = 0.0;
387 for(index =1; index < numberOfDaughters -1; index++) {
388 for(index2 = index+1; index2 < numberOfDaughters; index2++) {
389 if (rd[index] < rd[index2]){
390 temp = rd[index];
391 rd[index] = rd[index2];
392 rd[index2] = temp;
393 }
394 }
395 }
396
397 //calcurate virtual mass
398 tmas = parentmass - sumofdaughtermass;
399 temp = sumofdaughtermass;
400 for(index =0; index < numberOfDaughters; index++) {
401 sm[index] = rd[index]*tmas + temp;
402 temp -= daughtermass[index];
403 if (GetVerboseLevel()>1) {
404 G4cout << index << " rundom number:" << rd[index];
405 G4cout << " virtual mass:" << sm[index]/GeV << "[GeV/c/c]" <<G4endl;
406 }
407 }
408 delete [] rd;
409
410 //Calculate daughter momentum
411 weight = 1.0;
412 index =numberOfDaughters-1;
413 daughtermomentum[index]= Pmx( sm[index-1],daughtermass[index-1], sm[index]);
414#ifdef G4VERBOSE
415 if (GetVerboseLevel()>1) {
416 G4cout << " daughter " << index << ":" << *daughters_name[index];
417 G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
418 }
419#endif
420 for(index =numberOfDaughters-2; index>=0; index--) {
421 // calculate
422 daughtermomentum[index]= Pmx( sm[index],daughtermass[index], sm[index +1]);
423 if(daughtermomentum[index] < 0.0) {
424 // !!! illegal momentum !!!
425#ifdef G4VERBOSE
426 if (GetVerboseLevel()>0) {
427 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
428 G4cout << " can not calculate daughter momentum " <<G4endl;
429 G4cout << " parent:" << *parent_name;
430 G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
431 G4cout << " daughter " << index << ":" << *daughters_name[index];
432 G4cout << " mass:" << daughtermass[index]/GeV << "[GeV/c/c]" ;
433 G4cout << " mass:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
434 }
435#endif
436 delete [] sm;
437 delete [] daughtermass;
438 delete [] daughtermomentum;
439
440 G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt",
441 "PART112", JustWarning,
442 "Can not create decay products: sum of daughter mass is larger than parent mass");
443
444 return 0; // Error detection
445
446 } else {
447 // calculate weight of this events
448 weight *= daughtermomentum[index]/sm[index];
449#ifdef G4VERBOSE
450 if (GetVerboseLevel()>1) {
451 G4cout << " daughter " << index << ":" << *daughters_name[index];
452 G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
453 }
454#endif
455 }
456 }
457
458#ifdef G4VERBOSE
459 if (GetVerboseLevel()>1) {
460 G4cout << " weight: " << weight <<G4endl;
461 }
462#endif
463
464 // exit if number of Try exceeds 100
465 if (numberOfTry++ >100) {
466#ifdef G4VERBOSE
467 if (GetVerboseLevel()>0) {
468 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ";
469 G4cout << " can not determine Decay Kinematics " << G4endl;
470 G4cout << "parent : " << *parent_name << G4endl;
471 G4cout << "daughters : ";
472 for(index=0; index<numberOfDaughters; index++) {
473 G4cout << *daughters_name[index] << " , ";
474 }
475 G4cout << G4endl;
476 }
477#endif
478 G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ",
479 "PART113", JustWarning,
480 " Cannot decay : Decay Kinematics cannot be calculated ");
481
482 delete [] sm;
483 delete [] daughtermass;
484 delete [] daughtermomentum;
485 return 0; // Error detection
486 }
487 } while ( weight > G4UniformRand());
488
489#ifdef G4VERBOSE
490 if (GetVerboseLevel()>1) {
491 G4cout << "Start calulation of daughters momentum vector "<<G4endl;
492 }
493#endif
494
495 G4double costheta, sintheta, phi;
496 G4double beta;
497 daughterparticle = new G4DynamicParticle*[numberOfDaughters];
498
499 index = numberOfDaughters -2;
500 costheta = 2.*G4UniformRand()-1.0;
501 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
502 phi = twopi*G4UniformRand()*rad;
503 direction.setZ(costheta);
504 direction.setY(sintheta*std::sin(phi));
505 direction.setX(sintheta*std::cos(phi));
506 daughterparticle[index] = new G4DynamicParticle( daughters[index], direction*daughtermomentum[index] );
507 daughterparticle[index+1] = new G4DynamicParticle( daughters[index+1], direction*(-1.0*daughtermomentum[index]) );
508
509 for (index = numberOfDaughters -3; index >= 0; index--) {
510 //calculate momentum direction
511 costheta = 2.*G4UniformRand()-1.0;
512 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
513 phi = twopi*G4UniformRand()*rad;
514 direction.setZ(costheta);
515 direction.setY(sintheta*std::sin(phi));
516 direction.setX(sintheta*std::cos(phi));
517
518 // boost already created particles
519 beta = daughtermomentum[index];
520 beta /= std::sqrt( daughtermomentum[index]*daughtermomentum[index] + sm[index+1]*sm[index+1] );
521 for (index2 = index+1; index2<numberOfDaughters; index2++) {
523 // make G4LorentzVector for secondaries
524 p4 = daughterparticle[index2]->Get4Momentum();
525
526 // boost secondaries to new frame
527 p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
528
529 // change energy/momentum
530 daughterparticle[index2]->Set4Momentum(p4);
531 }
532 //create daughter G4DynamicParticle
533 daughterparticle[index]= new G4DynamicParticle( daughters[index], direction*(-1.0*daughtermomentum[index]));
534 }
535
536 //create G4Decayproducts
537 G4DynamicParticle *parentparticle;
538 direction.setX(1.0); direction.setY(0.0); direction.setZ(0.0);
539 parentparticle = new G4DynamicParticle( parent, direction, 0.0);
540 products = new G4DecayProducts(*parentparticle);
541 delete parentparticle;
542 for (index = 0; index<numberOfDaughters; index++) {
543 products->PushProducts(daughterparticle[index]);
544 }
545
546#ifdef G4VERBOSE
547 if (GetVerboseLevel()>1) {
548 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
549 G4cout << " create decay products in rest frame " << G4endl;
550 products->DumpInfo();
551 }
552#endif
553
554 delete [] daughterparticle;
555 delete [] daughtermomentum;
556 delete [] daughtermass;
557 delete [] sm;
558
559 return products;
560}
561
562
564{
565 // calcurate momentum of daughter particles in two-body decay
566 G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e);
567 if (ppp>0) return std::sqrt(ppp);
568 else return -1.;
569}
570
571
572
@ JustWarning
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cerr
G4DLLIMPORT std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:53
double z() const
double x() const
void setY(double)
double y() const
void setZ(double)
double mag() const
void setX(double)
HepLorentzVector & boost(double, double, double)
void DumpInfo() const
G4int PushProducts(G4DynamicParticle *aParticle)
G4LorentzVector Get4Momentum() const
void Set4Momentum(const G4LorentzVector &momentum)
const G4String & GetParticleName() const
static G4double Pmx(G4double e, G4double p1, G4double p2)
virtual G4DecayProducts * DecayIt(G4double)
G4String * parent_name
G4String ** daughters_name
G4int GetVerboseLevel() const
G4double * daughters_mass
G4ParticleDefinition * parent
G4ParticleDefinition ** daughters
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41