Geant4  10.01.p02
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 82335 2014-06-16 10:00:42Z 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, width
162  G4double daughtermass[2], daughterwidth[2];
163  G4double daughtermomentum;
164  daughtermass[0] = G4MT_daughters_mass[0];
165  daughtermass[1] = G4MT_daughters_mass[1];
166  daughterwidth[0] = G4MT_daughters_width[0];
167  daughterwidth[1] = G4MT_daughters_width[1];
168 
169  //create parent G4DynamicParticle at rest
170  G4ThreeVector dummy;
171  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
172  //create G4Decayproducts
173  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
174  delete parentparticle;
175 
176  G4bool withWidth = (daughterwidth[0]>1.0e-3*daughtermass[0]) || (daughterwidth[1]>1.0e-3*daughtermass[1]);
177  if (withWidth) {
178  G4double sumofdaughterwidthsq = daughterwidth[0]*daughterwidth[0]+daughterwidth[1]*daughterwidth[1];
179  // check parent mass and daughter mass
180  G4double maxDev = (parentmass - daughtermass[0] - daughtermass[1] )/std::sqrt(sumofdaughterwidthsq);
181  if (maxDev < -0.9*rangeMass ){
182 #ifdef G4VERBOSE
183  if (GetVerboseLevel()>0) {
184  G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt "
185  << "sum of daughter mass is larger than parent mass" << G4endl;
186  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass/GeV << G4endl;
187  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
188  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
189  }
190 #endif
191  G4Exception("G4PhaseSpaceDecayChannel::TwoBodyDecayIt",
192  "PART112", JustWarning,
193  "Can not create decay products: sum of daughter mass is larger than parent mass");
194  return products;
195  }
196  G4double dm1= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
197  G4double dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
198  while (dm1+dm2>parentmass){
199  dm1= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
200  dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
201  }
202  daughtermass[0] = dm1;
203  daughtermass[1] = dm2;
204  } else {
205  if (parentmass < daughtermass[0] + daughtermass[1] ){
206 #ifdef G4VERBOSE
207  if (GetVerboseLevel()>0) {
208  G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt "
209  << "sum of daughter mass is larger than parent mass" << G4endl;
210  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass/GeV << G4endl;
211  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
212  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
213  }
214 #endif
215  G4Exception("G4PhaseSpaceDecayChannel::TwoBodyDecayIt",
216  "PART112", JustWarning,
217  "Can not create decay products: sum of daughter mass is larger than parent mass");
218  return products;
219  }
220  }
221 
222  //calculate daughter momentum
223  daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
224 
225  G4double costheta = 2.*G4UniformRand()-1.0;
226  G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
227  G4double phi = twopi*G4UniformRand()*rad;
228  G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
229 
230  //create daughter G4DynamicParticle
231  G4double Ekin = std::sqrt(daughtermomentum*daughtermomentum + daughtermass[0]*daughtermass[0]) - daughtermass[0];
232  G4DynamicParticle * daughterparticle = new G4DynamicParticle( G4MT_daughters[0], direction, Ekin, daughtermass[0]);
233  products->PushProducts(daughterparticle);
234  Ekin = std::sqrt(daughtermomentum*daughtermomentum + daughtermass[1]*daughtermass[1]) - daughtermass[1];
235  daughterparticle = new G4DynamicParticle( G4MT_daughters[1], -1.0*direction, Ekin, daughtermass[1]);
236  products->PushProducts(daughterparticle);
237 
238 #ifdef G4VERBOSE
239  if (GetVerboseLevel()>1) {
240  G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt ";
241  G4cout << " create decay products in rest frame " <<G4endl;
242  products->DumpInfo();
243  }
244 #endif
245  return products;
246 }
247 
248 G4DecayProducts *G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()
249 // algorism of this code is originally written in GDECA3 of GEANT3
250 {
251 #ifdef G4VERBOSE
252  if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()"<<G4endl;
253 #endif
254  // parent mass
255  G4double parentmass = current_parent_mass;
256  //daughters'mass
257  G4double daughtermass[3], daughterwidth[3];
258  G4double sumofdaughtermass = 0.0;
259  G4double sumofdaughterwidthsq = 0.0;
260  G4bool withWidth = false;
261  for (G4int index=0; index<3; index++){
262  daughtermass[index] = G4MT_daughters_mass[index];
263  sumofdaughtermass += daughtermass[index];
264  daughterwidth[index] = G4MT_daughters_width[index];
265  sumofdaughterwidthsq += daughterwidth[index]*daughterwidth[index];
266  withWidth = withWidth ||(daughterwidth[index]>1.0e-3*daughtermass[index]);
267  }
268 
269  //create parent G4DynamicParticle at rest
270  G4ThreeVector dummy;
271  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
272 
273 
274  //create G4Decayproducts
275  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
276  delete parentparticle;
277 
278  if (withWidth){
279  G4double maxDev = (parentmass - sumofdaughtermass )/std::sqrt(sumofdaughterwidthsq) ;
280  if (maxDev < -0.9*rangeMass ){
281 #ifdef G4VERBOSE
282  if (GetVerboseLevel()>0) {
283  G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt "
284  << "sum of daughter mass is larger than parent mass" << G4endl;
285  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass/GeV << G4endl;
286  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
287  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
288  G4cout << "daughter 3:" << G4MT_daughters[2]->GetParticleName() << " " << daughtermass[2]/GeV << G4endl;
289  }
290 #endif
291  G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
292  "PART112", JustWarning,
293  "Can not create decay products: sum of daughter mass is larger than parent mass");
294  return products;
295  }
296  G4double dm1= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
297  G4double dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
298  G4double dm3= DynamicalMass(daughtermass[2],daughterwidth[2], maxDev);
299  while (dm1+dm2+dm3>parentmass){
300  dm1= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
301  dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
302  dm3= DynamicalMass(daughtermass[2],daughterwidth[2], maxDev);
303  }
304  daughtermass[0] = dm1;
305  daughtermass[1] = dm2;
306  daughtermass[2] = dm3;
307 
308  } else {
309  if (sumofdaughtermass >parentmass) {
310 #ifdef G4VERBOSE
311  if (GetVerboseLevel()>0) {
312  G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt "
313  << "sum of daughter mass is larger than parent mass" << G4endl;
314  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass/GeV << G4endl;
315  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
316  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
317  G4cout << "daughter 3:" << G4MT_daughters[2]->GetParticleName() << " " << daughtermass[2]/GeV << G4endl;
318  }
319 #endif
320  G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
321  "PART112", JustWarning,
322  "Can not create decay products: sum of daughter mass is larger than parent mass");
323  return products;
324  }
325  }
326 
327 
328 
329  //calculate daughter momentum
330  // Generate two
331  G4double rd1, rd2, rd;
332  G4double daughtermomentum[3];
333  G4double momentummax=0.0, momentumsum = 0.0;
334  G4double energy;
335 
336  do {
337  rd1 = G4UniformRand();
338  rd2 = G4UniformRand();
339  if (rd2 > rd1) {
340  rd = rd1;
341  rd1 = rd2;
342  rd2 = rd;
343  }
344  momentummax = 0.0;
345  momentumsum = 0.0;
346  // daughter 0
347  energy = rd2*(parentmass - sumofdaughtermass);
348  daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
349  if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0];
350  momentumsum += daughtermomentum[0];
351  // daughter 1
352  energy = (1.-rd1)*(parentmass - sumofdaughtermass);
353  daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
354  if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1];
355  momentumsum += daughtermomentum[1];
356  // daughter 2
357  energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
358  daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
359  if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2];
360  momentumsum += daughtermomentum[2];
361  } while (momentummax > momentumsum - momentummax );
362  // output message
363 #ifdef G4VERBOSE
364  if (GetVerboseLevel()>1) {
365  G4cout << " daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl;
366  G4cout << " daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl;
367  G4cout << " daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl;
368  G4cout << " momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl;
369  }
370 #endif
371 
372  //create daughter G4DynamicParticle
373  G4double costheta, sintheta, phi, sinphi, cosphi;
374  G4double costhetan, sinthetan, phin, sinphin, cosphin;
375  costheta = 2.*G4UniformRand()-1.0;
376  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
377  phi = twopi*G4UniformRand()*rad;
378  sinphi = std::sin(phi);
379  cosphi = std::cos(phi);
380 
381  G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
382  G4double Ekin = std::sqrt(daughtermomentum[0]*daughtermomentum[0] + daughtermass[0]*daughtermass[0]) - daughtermass[0];
383  G4DynamicParticle * daughterparticle = new G4DynamicParticle( G4MT_daughters[0],
384  direction0,
385  Ekin, daughtermass[0]);
386  products->PushProducts(daughterparticle);
387 
388  costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
389  sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
390  phin = twopi*G4UniformRand()*rad;
391  sinphin = std::sin(phin);
392  cosphin = std::cos(phin);
393  G4ThreeVector direction2;
394  direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
395  direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
396  direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
397  G4ThreeVector pmom = daughtermomentum[2]*direction2/direction2.mag();
398  Ekin = std::sqrt(pmom.mag2() + daughtermass[2]*daughtermass[2]) - daughtermass[2];
399  daughterparticle = new G4DynamicParticle( G4MT_daughters[2],
400  pmom/pmom.mag(),
401  Ekin, daughtermass[2]);
402  products->PushProducts(daughterparticle);
403 
404  pmom = (direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0);
405  Ekin = std::sqrt(pmom.mag2() + daughtermass[1]*daughtermass[1]) - daughtermass[1];
406  daughterparticle = new G4DynamicParticle(
407  G4MT_daughters[1],
408  pmom/pmom.mag(),
409  Ekin, daughtermass[1]);
410  products->PushProducts(daughterparticle);
411 
412 #ifdef G4VERBOSE
413  if (GetVerboseLevel()>1) {
414  G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt ";
415  G4cout << " create decay products in rest frame " <<G4endl;
416  products->DumpInfo();
417  }
418 #endif
419  return products;
420 }
421 
422 G4DecayProducts *G4PhaseSpaceDecayChannel::ManyBodyDecayIt()
423 // algorithm of this code is originally written in FORTRAN by M.Asai
424 //******************************************************************
425 // NBODY
426 // N-body phase space Monte-Carlo generator
427 // Makoto Asai
428 // Hiroshima Institute of Technology
429 // (asai@kekvax.kek.jp)
430 // Revised release : 19/Apr/1995
431 //
432 {
433  G4int index, index2;
434 
435  //return value
436  G4DecayProducts *products;
437 
438 #ifdef G4VERBOSE
439  if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()"<<G4endl;
440 #endif
441  // parent mass
442  G4double parentmass = current_parent_mass;
443  //daughters'mass
444  G4double *daughtermass = new G4double[numberOfDaughters];
445  G4double sumofdaughtermass = 0.0;
446  for (index=0; index<numberOfDaughters; index++){
447  daughtermass[index] = G4MT_daughters_mass[index];
448  sumofdaughtermass += daughtermass[index];
449  }
450 
451  //Calculate daughter momentum
452  G4double *daughtermomentum = new G4double[numberOfDaughters];
453  G4ThreeVector direction;
454  G4DynamicParticle **daughterparticle;
455  G4double *sm = new G4double[numberOfDaughters];
456  G4double tmas;
457  G4double weight = 1.0;
458  G4int numberOfTry = 0;
459 
460  do {
461  //Generate rundom number in descending order
462  G4double temp;
463  G4double *rd = new G4double[numberOfDaughters];
464  rd[0] = 1.0;
465  for(index =1; index < numberOfDaughters -1; index++) rd[index] = G4UniformRand();
466  rd[ numberOfDaughters -1] = 0.0;
467  for(index =1; index < numberOfDaughters -1; index++) {
468  for(index2 = index+1; index2 < numberOfDaughters; index2++) {
469  if (rd[index] < rd[index2]){
470  temp = rd[index];
471  rd[index] = rd[index2];
472  rd[index2] = temp;
473  }
474  }
475  }
476 
477  //calcurate virtual mass
478  tmas = parentmass - sumofdaughtermass;
479  temp = sumofdaughtermass;
480  for(index =0; index < numberOfDaughters; index++) {
481  sm[index] = rd[index]*tmas + temp;
482  temp -= daughtermass[index];
483  if (GetVerboseLevel()>1) {
484  G4cout << index << " rundom number:" << rd[index];
485  G4cout << " virtual mass:" << sm[index]/GeV << "[GeV/c/c]" <<G4endl;
486  }
487  }
488  delete [] rd;
489 
490  //Calculate daughter momentum
491  weight = 1.0;
492  index =numberOfDaughters-1;
493  daughtermomentum[index]= Pmx( sm[index-1],daughtermass[index-1], sm[index]);
494 #ifdef G4VERBOSE
495  if (GetVerboseLevel()>1) {
496  G4cout << " daughter " << index << ":" << *daughters_name[index];
497  G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
498  }
499 #endif
500  for(index =numberOfDaughters-2; index>=0; index--) {
501  // calculate
502  daughtermomentum[index]= Pmx( sm[index],daughtermass[index], sm[index +1]);
503  if(daughtermomentum[index] < 0.0) {
504  // !!! illegal momentum !!!
505 #ifdef G4VERBOSE
506  if (GetVerboseLevel()>0) {
507  G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
508  G4cout << " can not calculate daughter momentum " <<G4endl;
509  G4cout << " parent:" << *parent_name;
510  G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
511  G4cout << " daughter " << index << ":" << *daughters_name[index];
512  G4cout << " mass:" << daughtermass[index]/GeV << "[GeV/c/c]" ;
513  G4cout << " mass:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
514  }
515 #endif
516  delete [] sm;
517  delete [] daughtermass;
518  delete [] daughtermomentum;
519 
520  G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt",
521  "PART112", JustWarning,
522  "Can not create decay products: sum of daughter mass is larger than parent mass");
523 
524  return 0; // Error detection
525 
526  } else {
527  // calculate weight of this events
528  weight *= daughtermomentum[index]/sm[index];
529 #ifdef G4VERBOSE
530  if (GetVerboseLevel()>1) {
531  G4cout << " daughter " << index << ":" << *daughters_name[index];
532  G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
533  }
534 #endif
535  }
536  }
537 
538 #ifdef G4VERBOSE
539  if (GetVerboseLevel()>1) {
540  G4cout << " weight: " << weight <<G4endl;
541  }
542 #endif
543 
544  // exit if number of Try exceeds 100
545  if (numberOfTry++ >100) {
546 #ifdef G4VERBOSE
547  if (GetVerboseLevel()>0) {
548  G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ";
549  G4cout << " can not determine Decay Kinematics " << G4endl;
550  G4cout << "parent : " << *parent_name << G4endl;
551  G4cout << "daughters : ";
552  for(index=0; index<numberOfDaughters; index++) {
553  G4cout << *daughters_name[index] << " , ";
554  }
555  G4cout << G4endl;
556  }
557 #endif
558  G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ",
559  "PART113", JustWarning,
560  " Cannot decay : Decay Kinematics cannot be calculated ");
561 
562  delete [] sm;
563  delete [] daughtermass;
564  delete [] daughtermomentum;
565  return 0; // Error detection
566  }
567  } while ( weight > G4UniformRand());
568 
569 #ifdef G4VERBOSE
570  if (GetVerboseLevel()>1) {
571  G4cout << "Start calulation of daughters momentum vector "<<G4endl;
572  }
573 #endif
574 
575  G4double costheta, sintheta, phi;
576  G4double beta;
577  daughterparticle = new G4DynamicParticle*[numberOfDaughters];
578 
579  index = numberOfDaughters -2;
580  costheta = 2.*G4UniformRand()-1.0;
581  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
582  phi = twopi*G4UniformRand()*rad;
583  direction.setZ(costheta);
584  direction.setY(sintheta*std::sin(phi));
585  direction.setX(sintheta*std::cos(phi));
586  daughterparticle[index] = new G4DynamicParticle( G4MT_daughters[index], direction*daughtermomentum[index] );
587  daughterparticle[index+1] = new G4DynamicParticle( G4MT_daughters[index+1], direction*(-1.0*daughtermomentum[index]) );
588 
589  for (index = numberOfDaughters -3; index >= 0; index--) {
590  //calculate momentum direction
591  costheta = 2.*G4UniformRand()-1.0;
592  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
593  phi = twopi*G4UniformRand()*rad;
594  direction.setZ(costheta);
595  direction.setY(sintheta*std::sin(phi));
596  direction.setX(sintheta*std::cos(phi));
597 
598  // boost already created particles
599  beta = daughtermomentum[index];
600  beta /= std::sqrt( daughtermomentum[index]*daughtermomentum[index] + sm[index+1]*sm[index+1] );
601  for (index2 = index+1; index2<numberOfDaughters; index2++) {
602  G4LorentzVector p4;
603  // make G4LorentzVector for secondaries
604  p4 = daughterparticle[index2]->Get4Momentum();
605 
606  // boost secondaries to new frame
607  p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
608 
609  // change energy/momentum
610  daughterparticle[index2]->Set4Momentum(p4);
611  }
612  //create daughter G4DynamicParticle
613  daughterparticle[index]= new G4DynamicParticle( G4MT_daughters[index], direction*(-1.0*daughtermomentum[index]));
614  }
615 
616  //create G4Decayproducts
617  G4DynamicParticle *parentparticle;
618  direction.setX(1.0); direction.setY(0.0); direction.setZ(0.0);
619  parentparticle = new G4DynamicParticle( G4MT_parent, direction, 0.0);
620  products = new G4DecayProducts(*parentparticle);
621  delete parentparticle;
622  for (index = 0; index<numberOfDaughters; index++) {
623  products->PushProducts(daughterparticle[index]);
624  }
625 
626 #ifdef G4VERBOSE
627  if (GetVerboseLevel()>1) {
628  G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
629  G4cout << " create decay products in rest frame " << G4endl;
630  products->DumpInfo();
631  }
632 #endif
633 
634  delete [] daughterparticle;
635  delete [] daughtermomentum;
636  delete [] daughtermass;
637  delete [] sm;
638 
639  return products;
640 }
641 
642 
643 G4double G4PhaseSpaceDecayChannel::Pmx(G4double e, G4double p1, G4double p2)
644 {
645  // calcurate momentum of daughter particles in two-body decay
646  G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e);
647  if (ppp>0) return std::sqrt(ppp);
648  else return -1.;
649 }
650 
651 
652 
G4PhaseSpaceDecayChannel::DecayIt
virtual G4DecayProducts * DecayIt(G4double)
Definition: G4PhaseSpaceDecayChannel.cc:80
G4LorentzVector.hh
G4VDecayChannel::G4MT_daughters_mass
G4double * G4MT_daughters_mass
Definition: G4VDecayChannel.hh:245
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:242
G4VDecayChannel::G4MT_daughters
G4ParticleDefinition ** G4MT_daughters
Definition: G4VDecayChannel.hh:243
G4VDecayChannel::G4MT_parent_mass
G4double G4MT_parent_mass
Definition: G4VDecayChannel.hh:244
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:89
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:93
G4cout
G4GLOB_DLL std::ostream G4cout
G4VDecayChannel.hh
G4bool
bool G4bool
Definition: G4Types.hh:79
G4DecayProducts::DumpInfo
void DumpInfo() const
Definition: G4DecayProducts.cc:311
p2
const G4double p2
Definition: G4PreCompoundAlpha.cc:50
p1
const G4double p1
Definition: G4PreCompoundAlpha.cc:49
G4VDecayChannel::parent_name
G4String * parent_name
Definition: G4VDecayChannel.hh:220
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:376
G4VDecayChannel::GetVerboseLevel
G4int GetVerboseLevel() const
Definition: G4VDecayChannel.hh:351
G4VDecayChannel::G4MT_daughters_width
G4double * G4MT_daughters_width
Definition: G4VDecayChannel.hh:246
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:157
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VDecayChannel::DynamicalMass
G4double DynamicalMass(G4double massPDG, G4double width, G4double maxDev=+1.) const
Definition: G4VDecayChannel.cc:556
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4VDecayChannel::daughters_name
G4String ** daughters_name
Definition: G4VDecayChannel.hh:222
G4PhaseSpaceDecayChannel::Pmx
static G4double Pmx(G4double e, G4double p1, G4double p2)
Definition: G4PhaseSpaceDecayChannel.cc:643
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
G4LorentzRotation.hh