Geant4  10.00.p03
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: G4PhaseSpaceDecayChannel.cc 74261 2013-10-02 14:36:19Z gcosmo $
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 
39 #include "G4ParticleDefinition.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4DecayProducts.hh"
43 #include "G4VDecayChannel.hh"
44 #include "G4PhaseSpaceDecayChannel.hh"
45 #include "Randomize.hh"
46 #include "G4LorentzVector.hh"
47 #include "G4LorentzRotation.hh"
48 
49 G4ThreadLocal G4double G4PhaseSpaceDecayChannel::current_parent_mass = 0.0;
50 
51 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(G4int Verbose)
52  :G4VDecayChannel("Phase Space", Verbose)
53 {
54 
55 }
56 
57 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(
58  const G4String& theParentName,
59  G4double theBR,
60  G4int theNumberOfDaughters,
61  const G4String& theDaughterName1,
62  const G4String& theDaughterName2,
63  const G4String& theDaughterName3,
64  const G4String& theDaughterName4 )
65  :G4VDecayChannel("Phase Space",
66  theParentName,theBR,
67  theNumberOfDaughters,
68  theDaughterName1,
69  theDaughterName2,
70  theDaughterName3,
71  theDaughterName4)
72 {
73 
74 }
75 
76 G4PhaseSpaceDecayChannel::~G4PhaseSpaceDecayChannel()
77 {
78 }
79 
80 G4DecayProducts *G4PhaseSpaceDecayChannel::DecayIt(G4double parentMass)
81 {
82 #ifdef G4VERBOSE
83  if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
84 #endif
85 
86  G4DecayProducts * products = 0;
87 
88  if (G4MT_parent == 0) FillParent();
89  if (G4MT_daughters == 0) FillDaughters();
90 
91  if (parentMass >0.0) current_parent_mass = parentMass;
92  else current_parent_mass = G4MT_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 
126 G4DecayProducts *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( G4MT_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( G4MT_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 
153 G4DecayProducts *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] = G4MT_daughters_mass[0];
165  daughtermass[1] = G4MT_daughters_mass[1];
166 
167  //create parent G4DynamicParticle at rest
168  G4ThreeVector dummy;
169  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_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  G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt "
180  << "sum of daughter mass is larger than parent mass" << G4endl;
181  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass/GeV << G4endl;
182  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
183  G4cout << "daughter 2:" << G4MT_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( G4MT_daughters[0], direction*daughtermomentum);
199  products->PushProducts(daughterparticle);
200  daughterparticle = new G4DynamicParticle( G4MT_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 
213 G4DecayProducts *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] = G4MT_daughters_mass[index];
226  sumofdaughtermass += daughtermass[index];
227  }
228 
229  //create parent G4DynamicParticle at rest
230  G4ThreeVector dummy;
231  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_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  G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt "
242  << "sum of daughter mass is larger than parent mass" << G4endl;
243  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass/GeV << G4endl;
244  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
245  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
246  G4cout << "daughter 3:" << G4MT_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( G4MT_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( G4MT_daughters[2], direction2*(daughtermomentum[2]/direction2.mag()));
323  products->PushProducts(daughterparticle);
324 
325  daughterparticle =
326  new G4DynamicParticle(
327  G4MT_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 
342 G4DecayProducts *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
349 // (asai@kekvax.kek.jp)
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] = G4MT_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;
375  G4double *sm = new G4double[numberOfDaughters];
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;
383  G4double *rd = new G4double[numberOfDaughters];
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( G4MT_daughters[index], direction*daughtermomentum[index] );
507  daughterparticle[index+1] = new G4DynamicParticle( G4MT_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++) {
522  G4LorentzVector p4;
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( G4MT_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( G4MT_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 
563 G4double G4PhaseSpaceDecayChannel::Pmx(G4double e, G4double p1, G4double p2)
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 
G4PhaseSpaceDecayChannel::DecayIt
virtual G4DecayProducts * DecayIt(G4double)
Definition: G4PhaseSpaceDecayChannel.cc:80
G4LorentzVector.hh
G4VDecayChannel::G4MT_daughters_mass
G4double * G4MT_daughters_mass
Definition: G4VDecayChannel.hh:233
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4DecayProducts::PushProducts
G4int PushProducts(G4DynamicParticle *aParticle)
Definition: G4DecayProducts.cc:153
G4VDecayChannel::G4MT_parent
G4ParticleDefinition * G4MT_parent
Definition: G4VDecayChannel.hh:230
G4VDecayChannel::G4MT_daughters
G4ParticleDefinition ** G4MT_daughters
Definition: G4VDecayChannel.hh:231
G4VDecayChannel::G4MT_parent_mass
G4double G4MT_parent_mass
Definition: G4VDecayChannel.hh:232
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:52
JustWarning
Definition: G4ExceptionSeverity.hh:64
G4int
int G4int
Definition: G4Types.hh:78
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:159
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4cout
G4GLOB_DLL std::ostream G4cout
G4VDecayChannel.hh
G4DecayProducts::DumpInfo
void DumpInfo() const
Definition: G4DecayProducts.cc:311
G4VDecayChannel::parent_name
G4String * parent_name
Definition: G4VDecayChannel.hh:211
GeV
static const double GeV
Definition: G4SIunits.hh:196
G4ParticleDefinition.hh
Randomize.hh
G4DynamicParticle::Get4Momentum
G4LorentzVector Get4Momentum() const
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4PhaseSpaceDecayChannel.hh
G4PhysicalConstants.hh
rad
static const double rad
Definition: G4SIunits.hh:130
G4DecayProducts.hh
G4PhaseSpaceDecayChannel::current_parent_mass
static G4ThreadLocal G4double current_parent_mass
Definition: G4PhaseSpaceDecayChannel.hh:69
G4DynamicParticle::Set4Momentum
void Set4Momentum(const G4LorentzVector &momentum)
Definition: G4DynamicParticle.cc:353
G4VDecayChannel::GetVerboseLevel
G4int GetVerboseLevel() const
Definition: G4VDecayChannel.hh:338
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:156
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4VDecayChannel::daughters_name
G4String ** daughters_name
Definition: G4VDecayChannel.hh:213
G4PhaseSpaceDecayChannel::Pmx
static G4double Pmx(G4double e, G4double p1, G4double p2)
Definition: G4PhaseSpaceDecayChannel.cc:563
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
G4LorentzRotation.hh