ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PhaseSpaceDecayChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PhaseSpaceDecayChannel.cc
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 //
29 // ------------------------------------------------------------
30 // GEANT 4 class header file
31 //
32 // History: first implementation, based on object model of
33 // 27 July 1996 H.Kurashige
34 // 20 Oct 1996 H.Kurashige
35 // 30 May 1997 H.Kurashige
36 // ------------------------------------------------------------
37 
38 #include "G4ParticleDefinition.hh"
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4DecayProducts.hh"
42 #include "G4VDecayChannel.hh"
43 #include "G4PhaseSpaceDecayChannel.hh"
44 #include "Randomize.hh"
45 #include "G4LorentzVector.hh"
46 #include "G4LorentzRotation.hh"
47 
48 //G4ThreadLocal G4double G4PhaseSpaceDecayChannel::current_parent_mass = 0.0;
49 
50 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(G4int Verbose)
51  :G4VDecayChannel("Phase Space", Verbose),
52  useGivenDaughterMass(false)
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  useGivenDaughterMass(false)
73 {
74 
75 }
76 
77 G4PhaseSpaceDecayChannel::~G4PhaseSpaceDecayChannel()
78 {
79 }
80 
81 G4DecayProducts *G4PhaseSpaceDecayChannel::DecayIt(G4double parentMass)
82 {
83 #ifdef G4VERBOSE
84  if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
85 #endif
86 
87  G4DecayProducts * products = 0;
88 
89  CheckAndFillParent();
90  CheckAndFillDaughters();
91 
92  if (parentMass >0.0) current_parent_mass.Put( parentMass );
93  else current_parent_mass.Put(G4MT_parent_mass);
94 
95  switch (numberOfDaughters){
96  case 0:
97 #ifdef G4VERBOSE
98  if (GetVerboseLevel()>0) {
99  G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
100  G4cout << " daughters not defined " <<G4endl;
101  }
102 #endif
103  break;
104  case 1:
105  products = OneBodyDecayIt();
106  break;
107  case 2:
108  products = TwoBodyDecayIt();
109  break;
110  case 3:
111  products = ThreeBodyDecayIt();
112  break;
113  default:
114  products = ManyBodyDecayIt();
115  break;
116  }
117 #ifdef G4VERBOSE
118  if ((products == 0) && (GetVerboseLevel()>0)) {
119  G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
120  G4cout << *parent_name << " can not decay " << G4endl;
121  DumpInfo();
122  }
123 #endif
124  return products;
125 }
126 
127 G4DecayProducts *G4PhaseSpaceDecayChannel::OneBodyDecayIt()
128 {
129 #ifdef G4VERBOSE
130  if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt()"<<G4endl;
131 #endif
132  // parent mass
133  G4double parentmass = current_parent_mass.Get();
134 
135  //create parent G4DynamicParticle at rest
136  G4ThreeVector dummy;
137  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0, parentmass);
138  //create G4Decayproducts
139  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
140  delete parentparticle;
141 
142  //create daughter G4DynamicParticle at rest
143  G4DynamicParticle * daughterparticle = new G4DynamicParticle( G4MT_daughters[0], dummy, 0.0);
144  if (useGivenDaughterMass) daughterparticle->SetMass(givenDaughterMasses[0]);
145  products->PushProducts(daughterparticle);
146 
147 #ifdef G4VERBOSE
148  if (GetVerboseLevel()>1) {
149  G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt ";
150  G4cout << " create decay products in rest frame " <<G4endl;
151  products->DumpInfo();
152  }
153 #endif
154  return products;
155 }
156 
157 G4DecayProducts *G4PhaseSpaceDecayChannel::TwoBodyDecayIt()
158 {
159 #ifdef G4VERBOSE
160  if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt()"<<G4endl;
161 #endif
162  // parent mass
163  G4double parentmass = current_parent_mass.Get();
164 
165  //daughters'mass, width
166  G4double daughtermass[2], daughterwidth[2];
167  daughtermass[0] = G4MT_daughters_mass[0];
168  daughtermass[1] = G4MT_daughters_mass[1];
169  daughterwidth[0] = G4MT_daughters_width[0];
170  daughterwidth[1] = G4MT_daughters_width[1];
171 
172  //create parent G4DynamicParticle at rest
173  G4ThreeVector dummy;
174  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0, parentmass);
175  //create G4Decayproducts
176  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
177  delete parentparticle;
178 
179  if (!useGivenDaughterMass) {
180  G4bool withWidth = (daughterwidth[0]>1.0e-3*daughtermass[0])
181  || (daughterwidth[1]>1.0e-3*daughtermass[1]);
182  if (withWidth) {
183  G4double sumofdaughterwidthsq = daughterwidth[0]*daughterwidth[0]+daughterwidth[1]*daughterwidth[1];
184  // check parent mass and daughter mass
185  G4double maxDev = (parentmass - daughtermass[0] - daughtermass[1] )/std::sqrt(sumofdaughterwidthsq);
186  if (maxDev <= -1.0*rangeMass ){
187 #ifdef G4VERBOSE
188  if (GetVerboseLevel()>0) {
189  G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt "
190  << "sum of daughter mass is larger than parent mass" << G4endl;
191  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass.Get()/GeV << G4endl;
192  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
193  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
194  }
195 #endif
196  G4Exception("G4PhaseSpaceDecayChannel::TwoBodyDecayIt",
197  "PART112", JustWarning,
198  "Can not create decay products: sum of daughter mass is larger than parent mass");
199  return products;
200  }
201  G4double dm1=daughtermass[0];
202  if (daughterwidth[0] > 0.) dm1= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
203  G4double dm2= daughtermass[1];
204  if (daughterwidth[1] > 0.) dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
205  while (dm1+dm2>parentmass){ // Loop checking, 09.08.2015, K.Kurashige
206  dm1= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
207  dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
208  }
209  daughtermass[0] = dm1;
210  daughtermass[1] = dm2;
211  }
212  } else {
213  // use given daughter mass;
214  daughtermass[0] = givenDaughterMasses[0];
215  daughtermass[1] = givenDaughterMasses[1];
216  }
217  if (parentmass < daughtermass[0] + daughtermass[1] ){
218 #ifdef G4VERBOSE
219  if (GetVerboseLevel()>0) {
220  G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt "
221  << "sum of daughter mass is larger than parent mass" << G4endl;
222  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass.Get()/GeV << G4endl;
223  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
224  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
225  if (useGivenDaughterMass) {
226  G4cout << "Daughter Mass is given " << G4endl;
227  }
228  }
229 #endif
230  G4Exception("G4PhaseSpaceDecayChannel::TwoBodyDecayIt",
231  "PART112", JustWarning,
232  "Can not create decay products: sum of daughter mass is larger than parent mass");
233  return products;
234  }
235 
236  //calculate daughter momentum
237  G4double daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
238 
239  G4double costheta = 2.*G4UniformRand()-1.0;
240  G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
241  G4double phi = twopi*G4UniformRand()*rad;
242  G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
243 
244  //create daughter G4DynamicParticle
245  G4double Ekin = std::sqrt(daughtermomentum*daughtermomentum + daughtermass[0]*daughtermass[0]) - daughtermass[0];
246  G4DynamicParticle * daughterparticle = new G4DynamicParticle( G4MT_daughters[0], direction, Ekin, daughtermass[0]);
247  products->PushProducts(daughterparticle);
248  Ekin = std::sqrt(daughtermomentum*daughtermomentum + daughtermass[1]*daughtermass[1]) - daughtermass[1];
249  daughterparticle = new G4DynamicParticle( G4MT_daughters[1], -1.0*direction, Ekin, daughtermass[1]);
250  products->PushProducts(daughterparticle);
251 
252 #ifdef G4VERBOSE
253  if (GetVerboseLevel()>1) {
254  G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt ";
255  G4cout << " create decay products in rest frame " <<G4endl;
256  products->DumpInfo();
257  }
258 #endif
259  return products;
260 }
261 
262 G4DecayProducts *G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()
263 // algorism of this code is originally written in GDECA3 of GEANT3
264 {
265 #ifdef G4VERBOSE
266  if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()"<<G4endl;
267 #endif
268  // parent mass
269  G4double parentmass = current_parent_mass.Get();
270  //daughters'mass
271  G4double daughtermass[3], daughterwidth[3];
272  G4double sumofdaughtermass = 0.0;
273  G4double sumofdaughterwidthsq = 0.0;
274  G4bool withWidth = false;
275  for (G4int index=0; index<3; index++){
276  daughtermass[index] = G4MT_daughters_mass[index];
277  sumofdaughtermass += daughtermass[index];
278  daughterwidth[index] = G4MT_daughters_width[index];
279  sumofdaughterwidthsq += daughterwidth[index]*daughterwidth[index];
280  withWidth = withWidth ||(daughterwidth[index]>1.0e-3*daughtermass[index]);
281  }
282 
283  //create parent G4DynamicParticle at rest
284  G4ThreeVector dummy;
285  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0, parentmass);
286 
287 
288  //create G4Decayproducts
289  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
290  delete parentparticle;
291 
292  if (!useGivenDaughterMass) {
293  if (withWidth){
294  G4double maxDev = (parentmass - sumofdaughtermass )/std::sqrt(sumofdaughterwidthsq) ;
295  if (maxDev <= -1.0*rangeMass ){
296 #ifdef G4VERBOSE
297  if (GetVerboseLevel()>0) {
298  G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt "
299  << "sum of daughter mass is larger than parent mass" << G4endl;
300  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass.Get()/GeV << G4endl;
301  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
302  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
303  G4cout << "daughter 3:" << G4MT_daughters[2]->GetParticleName() << " " << daughtermass[2]/GeV << G4endl;
304  }
305 #endif
306  G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
307  "PART112", JustWarning,
308  "Can not create decay products: sum of daughter mass is larger than parent mass");
309  return products;
310  }
311  G4double dm1=daughtermass[0];
312  if (daughterwidth[0] > 0.) dm1= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
313  G4double dm2= daughtermass[1];
314  if (daughterwidth[1] > 0.) dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
315  G4double dm3= daughtermass[2];
316  if (daughterwidth[2] > 0.) dm3= DynamicalMass(daughtermass[2],daughterwidth[2], maxDev);
317  while (dm1+dm2+dm3>parentmass){ // Loop checking, 09.08.2015, K.Kurashige
318  dm1= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
319  dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
320  dm3= DynamicalMass(daughtermass[2],daughterwidth[2], maxDev);
321  }
322  daughtermass[0] = dm1;
323  daughtermass[1] = dm2;
324  daughtermass[2] = dm3;
325  sumofdaughtermass = dm1+dm2+dm3;
326  }
327  } else {
328  // use given daughter mass;
329  daughtermass[0] = givenDaughterMasses[0];
330  daughtermass[1] = givenDaughterMasses[1];
331  daughtermass[2] = givenDaughterMasses[2];
332  sumofdaughtermass = daughtermass[0] + daughtermass[1] + daughtermass[2];
333  }
334 
335  if (sumofdaughtermass >parentmass) {
336 #ifdef G4VERBOSE
337  if (GetVerboseLevel()>0) {
338  G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt "
339  << "sum of daughter mass is larger than parent mass" << G4endl;
340  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass.Get()/GeV << G4endl;
341  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
342  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
343  G4cout << "daughter 3:" << G4MT_daughters[2]->GetParticleName() << " " << daughtermass[2]/GeV << G4endl;
344  }
345  if (useGivenDaughterMass) {
346  G4cout << "Daughter Mass is given " << G4endl;
347  }
348 #endif
349  G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
350  "PART112", JustWarning,
351  "Can not create decay products: sum of daughter mass is larger than parent mass");
352  return products;
353  }
354 
355 
356 
357  //calculate daughter momentum
358  // Generate two
359  G4double rd1, rd2, rd;
360  G4double daughtermomentum[3];
361  G4double momentummax=0.0, momentumsum = 0.0;
362  G4double energy;
363  const size_t MAX_LOOP=10000;
364 
365  for (size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
366  rd1 = G4UniformRand();
367  rd2 = G4UniformRand();
368  if (rd2 > rd1) {
369  rd = rd1;
370  rd1 = rd2;
371  rd2 = rd;
372  }
373  momentummax = 0.0;
374  momentumsum = 0.0;
375  // daughter 0
376  energy = rd2*(parentmass - sumofdaughtermass);
377  daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
378  if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0];
379  momentumsum += daughtermomentum[0];
380  // daughter 1
381  energy = (1.-rd1)*(parentmass - sumofdaughtermass);
382  daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
383  if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1];
384  momentumsum += daughtermomentum[1];
385  // daughter 2
386  energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
387  daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
388  if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2];
389  momentumsum += daughtermomentum[2];
390  if (momentummax <= momentumsum - momentummax ) break;
391  }
392 
393  // output message
394 #ifdef G4VERBOSE
395  if (GetVerboseLevel()>1) {
396  G4cout << " daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl;
397  G4cout << " daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl;
398  G4cout << " daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl;
399  G4cout << " momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl;
400  }
401 #endif
402 
403  //create daughter G4DynamicParticle
404  G4double costheta, sintheta, phi, sinphi, cosphi;
405  G4double costhetan, sinthetan, phin, sinphin, cosphin;
406  costheta = 2.*G4UniformRand()-1.0;
407  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
408  phi = twopi*G4UniformRand()*rad;
409  sinphi = std::sin(phi);
410  cosphi = std::cos(phi);
411 
412  G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
413  G4double Ekin = std::sqrt(daughtermomentum[0]*daughtermomentum[0] + daughtermass[0]*daughtermass[0]) - daughtermass[0];
414  G4DynamicParticle * daughterparticle = new G4DynamicParticle( G4MT_daughters[0],
415  direction0,
416  Ekin, daughtermass[0]);
417  products->PushProducts(daughterparticle);
418 
419  costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
420  sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
421  phin = twopi*G4UniformRand()*rad;
422  sinphin = std::sin(phin);
423  cosphin = std::cos(phin);
424  G4ThreeVector direction2;
425  direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
426  direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
427  direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
428  G4ThreeVector pmom = daughtermomentum[2]*direction2/direction2.mag();
429  Ekin = std::sqrt(pmom.mag2() + daughtermass[2]*daughtermass[2]) - daughtermass[2];
430  daughterparticle = new G4DynamicParticle( G4MT_daughters[2],
431  pmom/pmom.mag(),
432  Ekin, daughtermass[2]);
433  products->PushProducts(daughterparticle);
434 
435  pmom = (direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0);
436  Ekin = std::sqrt(pmom.mag2() + daughtermass[1]*daughtermass[1]) - daughtermass[1];
437  daughterparticle = new G4DynamicParticle(
438  G4MT_daughters[1],
439  pmom/pmom.mag(),
440  Ekin, daughtermass[1]);
441  products->PushProducts(daughterparticle);
442 
443 #ifdef G4VERBOSE
444  if (GetVerboseLevel()>1) {
445  G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt ";
446  G4cout << " create decay products in rest frame " <<G4endl;
447  products->DumpInfo();
448  }
449 #endif
450  return products;
451 }
452 
453 G4DecayProducts *G4PhaseSpaceDecayChannel::ManyBodyDecayIt()
454 // algorithm of this code is originally written in FORTRAN by M.Asai
455 //******************************************************************
456 // NBODY
457 // N-body phase space Monte-Carlo generator
458 // Makoto Asai
459 // Hiroshima Institute of Technology
460 // (asai@kekvax.kek.jp)
461 // Revised release : 19/Apr/1995
462 //
463 {
464  G4int index, index2;
465 
466 #ifdef G4VERBOSE
467  if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()"<<G4endl;
468 #endif
469 
470 
471  // parent mass
472  G4double parentmass = current_parent_mass.Get();
473 
474  //parent particle
475  G4ThreeVector dummy;
476  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0, parentmass);
477 
478  //create G4Decayproducts
479  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
480  delete parentparticle;
481 
482  //daughters'mass
483  G4double *daughtermass = new G4double[numberOfDaughters];
484 
485  G4double sumofdaughtermass = 0.0;
486  for (index=0; index<numberOfDaughters; index++){
487  if (!useGivenDaughterMass) {
488  daughtermass[index] = G4MT_daughters_mass[index];
489  } else {
490  daughtermass[index] = givenDaughterMasses[index];
491  }
492  sumofdaughtermass += daughtermass[index];
493  }
494  if (sumofdaughtermass >parentmass) {
495 #ifdef G4VERBOSE
496  if (GetVerboseLevel()>0) {
497  G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt "
498  << "sum of daughter mass is larger than parent mass" << G4endl;
499  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass.Get()/GeV << G4endl;
500  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
501  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
502  G4cout << "daughter 3:" << G4MT_daughters[2]->GetParticleName() << " " << daughtermass[2]/GeV << G4endl;
503  G4cout << "daughter 4:" << G4MT_daughters[3]->GetParticleName() << " " << daughtermass[3]/GeV << G4endl;
504  }
505 #endif
506  G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt",
507  "PART112", JustWarning,
508  "Can not create decay products: sum of daughter mass is larger than parent mass");
509  delete [] daughtermass;
510  return products;
511  }
512 
513  //Calculate daughter momentum
514  G4double *daughtermomentum = new G4double[numberOfDaughters];
515  G4ThreeVector direction;
516  G4DynamicParticle **daughterparticle;
517  G4double *sm = new G4double[numberOfDaughters];
518  G4double tmas;
519  G4double weight = 1.0;
520  G4int numberOfTry = 0;
521  const size_t MAX_LOOP=10000;
522 
523  for (size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
524  //Generate rundom number in descending order
525  G4double temp;
526  G4double *rd = new G4double[numberOfDaughters];
527  rd[0] = 1.0;
528  for(index =1; index < numberOfDaughters -1; index++) rd[index] = G4UniformRand();
529  rd[ numberOfDaughters -1] = 0.0;
530  for(index =1; index < numberOfDaughters -1; index++) {
531  for(index2 = index+1; index2 < numberOfDaughters; index2++) {
532  if (rd[index] < rd[index2]){
533  temp = rd[index];
534  rd[index] = rd[index2];
535  rd[index2] = temp;
536  }
537  }
538  }
539 
540  //calcurate virtual mass
541  tmas = parentmass - sumofdaughtermass;
542  temp = sumofdaughtermass;
543  for(index =0; index < numberOfDaughters; index++) {
544  sm[index] = rd[index]*tmas + temp;
545  temp -= daughtermass[index];
546  if (GetVerboseLevel()>1) {
547  G4cout << index << " rundom number:" << rd[index];
548  G4cout << " virtual mass:" << sm[index]/GeV << "[GeV/c/c]" <<G4endl;
549  }
550  }
551  delete [] rd;
552 
553  //Calculate daughter momentum
554  weight = 1.0;
555  G4bool smOK=true;
556  for(index =0; index< numberOfDaughters-1 && smOK; index++) {
557  smOK = (sm[index]-daughtermass[index]- sm[index+1] >=0.);
558  }
559  if (!smOK) continue;
560 
561  index =numberOfDaughters-1;
562  daughtermomentum[index]= Pmx( sm[index-1],daughtermass[index-1], sm[index]);
563 #ifdef G4VERBOSE
564  if (GetVerboseLevel()>1) {
565  G4cout << " daughter " << index << ":" << *daughters_name[index];
566  G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
567  }
568 #endif
569  for(index =numberOfDaughters-2; index>=0; index--) {
570  // calculate
571  daughtermomentum[index]= Pmx( sm[index],daughtermass[index], sm[index +1]);
572  if(daughtermomentum[index] < 0.0) {
573  // !!! illegal momentum !!!
574 #ifdef G4VERBOSE
575  if (GetVerboseLevel()>0) {
576  G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
577  G4cout << " can not calculate daughter momentum " <<G4endl;
578  G4cout << " parent:" << *parent_name;
579  G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
580  G4cout << " daughter " << index << ":" << *daughters_name[index];
581  G4cout << " mass:" << daughtermass[index]/GeV << "[GeV/c/c]" ;
582  G4cout << " mass:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
583  if (useGivenDaughterMass) {
584  G4cout << "Daughter Mass is given " << G4endl;
585  }
586  }
587 #endif
588  delete [] sm;
589  delete [] daughtermass;
590  delete [] daughtermomentum;
591  delete products;
592  G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt",
593  "PART112", JustWarning,
594  "Can not create decay products: sum of daughter mass is larger than parent mass");
595 
596  return 0; // Error detection
597 
598  } else {
599  // calculate weight of this events
600  weight *= daughtermomentum[index]/sm[index];
601 #ifdef G4VERBOSE
602  if (GetVerboseLevel()>1) {
603  G4cout << " daughter " << index << ":" << *daughters_name[index];
604  G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
605  }
606 #endif
607  }
608  }
609 
610 #ifdef G4VERBOSE
611  if (GetVerboseLevel()>1) {
612  G4cout << " weight: " << weight <<G4endl;
613  }
614 #endif
615 
616  // exit if number of Try exceeds 100
617  if (numberOfTry++ >100) {
618 #ifdef G4VERBOSE
619  if (GetVerboseLevel()>0) {
620  G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ";
621  G4cout << " can not determine Decay Kinematics " << G4endl;
622  G4cout << "parent : " << *parent_name << G4endl;
623  G4cout << "daughters : ";
624  for(index=0; index<numberOfDaughters; index++) {
625  G4cout << *daughters_name[index] << " , ";
626  }
627  G4cout << G4endl;
628  }
629 #endif
630  G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ",
631  "PART113", JustWarning,
632  " Cannot decay : Decay Kinematics cannot be calculated ");
633 
634  delete [] sm;
635  delete [] daughtermass;
636  delete [] daughtermomentum;
637  delete products;
638  return 0; // Error detection
639  }
640  if ( weight < G4UniformRand()) break;
641  }
642 
643 #ifdef G4VERBOSE
644  if (GetVerboseLevel()>1) {
645  G4cout << "Start calulation of daughters momentum vector "<<G4endl;
646  }
647 #endif
648 
649  G4double costheta, sintheta, phi;
650  G4double beta;
651  daughterparticle = new G4DynamicParticle*[numberOfDaughters];
652 
653  index = numberOfDaughters -2;
654  costheta = 2.*G4UniformRand()-1.0;
655  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
656  phi = twopi*G4UniformRand()*rad;
657  direction.setZ(costheta);
658  direction.setY(sintheta*std::sin(phi));
659  direction.setX(sintheta*std::cos(phi));
660  daughterparticle[index] = new G4DynamicParticle( G4MT_daughters[index], direction*daughtermomentum[index] );
661  daughterparticle[index+1] = new G4DynamicParticle( G4MT_daughters[index+1], direction*(-1.0*daughtermomentum[index]) );
662 
663  for (index = numberOfDaughters -3; index >= 0; index--) {
664  //calculate momentum direction
665  costheta = 2.*G4UniformRand()-1.0;
666  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
667  phi = twopi*G4UniformRand()*rad;
668  direction.setZ(costheta);
669  direction.setY(sintheta*std::sin(phi));
670  direction.setX(sintheta*std::cos(phi));
671 
672  // boost already created particles
673  beta = daughtermomentum[index];
674  beta /= std::sqrt( daughtermomentum[index]*daughtermomentum[index] + sm[index+1]*sm[index+1] );
675  for (index2 = index+1; index2<numberOfDaughters; index2++) {
676  G4LorentzVector p4;
677  // make G4LorentzVector for secondaries
678  p4 = daughterparticle[index2]->Get4Momentum();
679 
680  // boost secondaries to new frame
681  p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
682 
683  // change energy/momentum
684  daughterparticle[index2]->Set4Momentum(p4);
685  }
686  //create daughter G4DynamicParticle
687  daughterparticle[index]= new G4DynamicParticle( G4MT_daughters[index], direction*(-1.0*daughtermomentum[index]));
688  }
689 
690  //add daughters to G4Decayproducts
691  for (index = 0; index<numberOfDaughters; index++) {
692  products->PushProducts(daughterparticle[index]);
693  }
694 
695 #ifdef G4VERBOSE
696  if (GetVerboseLevel()>1) {
697  G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
698  G4cout << " create decay products in rest frame " << G4endl;
699  products->DumpInfo();
700  }
701 #endif
702 
703  delete [] daughterparticle;
704  delete [] daughtermomentum;
705  delete [] daughtermass;
706  delete [] sm;
707 
708  return products;
709 }
710 
711 G4bool G4PhaseSpaceDecayChannel::SetDaughterMasses( G4double masses[])
712 {
713  for (G4int idx=0; idx<numberOfDaughters; idx++){
714  givenDaughterMasses[idx] = masses[idx];
715  }
716  useGivenDaughterMass = true;
717  return useGivenDaughterMass;
718 }
719 
720 G4bool G4PhaseSpaceDecayChannel::SampleDaughterMasses()
721 {
722  useGivenDaughterMass = false;
723  return useGivenDaughterMass;
724 }
725 
726 G4bool G4PhaseSpaceDecayChannel::IsOKWithParentMass(G4double parentMass) {
727  if (!useGivenDaughterMass) return G4VDecayChannel::IsOKWithParentMass(parentMass);
728 
729  CheckAndFillParent();
730  CheckAndFillDaughters();
731 
732  G4double sumOfDaughterMassMin=0.0;
733  for (G4int index=0; index < numberOfDaughters; index++) {
734  sumOfDaughterMassMin += givenDaughterMasses[index];
735  }
736  return (parentMass >= sumOfDaughterMassMin);
737 }
738 
739 G4double G4PhaseSpaceDecayChannel::Pmx(G4double e, G4double p1, G4double p2)
740 {
741  // calcurate momentum of daughter particles in two-body decay
742  G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e);
743  if (ppp>0) return std::sqrt(ppp);
744  else return -1.;
745 }
746 
747 
748