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