Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 98352 2016-07-08 08:21:00Z 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  useGivenDaughterMass(false)
54 {
55 
56 }
57 
58 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(
59  const G4String& theParentName,
60  G4double theBR,
61  G4int theNumberOfDaughters,
62  const G4String& theDaughterName1,
63  const G4String& theDaughterName2,
64  const G4String& theDaughterName3,
65  const G4String& theDaughterName4 )
66  :G4VDecayChannel("Phase Space",
67  theParentName,theBR,
68  theNumberOfDaughters,
69  theDaughterName1,
70  theDaughterName2,
71  theDaughterName3,
72  theDaughterName4),
73  useGivenDaughterMass(false)
74 {
75 
76 }
77 
78 G4PhaseSpaceDecayChannel::~G4PhaseSpaceDecayChannel()
79 {
80 }
81 
82 G4DecayProducts *G4PhaseSpaceDecayChannel::DecayIt(G4double parentMass)
83 {
84 #ifdef G4VERBOSE
85  if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
86 #endif
87 
88  G4DecayProducts * products = 0;
89 
90  CheckAndFillParent();
91  CheckAndFillDaughters();
92 
93  if (parentMass >0.0) current_parent_mass.Put( parentMass );
94  else current_parent_mass.Put(G4MT_parent_mass);
95 
96  switch (numberOfDaughters){
97  case 0:
98 #ifdef G4VERBOSE
99  if (GetVerboseLevel()>0) {
100  G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
101  G4cout << " daughters not defined " <<G4endl;
102  }
103 #endif
104  break;
105  case 1:
106  products = OneBodyDecayIt();
107  break;
108  case 2:
109  products = TwoBodyDecayIt();
110  break;
111  case 3:
112  products = ThreeBodyDecayIt();
113  break;
114  default:
115  products = ManyBodyDecayIt();
116  break;
117  }
118 #ifdef G4VERBOSE
119  if ((products == 0) && (GetVerboseLevel()>0)) {
120  G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
121  G4cout << *parent_name << " can not decay " << G4endl;
122  DumpInfo();
123  }
124 #endif
125  return products;
126 }
127 
128 G4DecayProducts *G4PhaseSpaceDecayChannel::OneBodyDecayIt()
129 {
130 #ifdef G4VERBOSE
131  if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt()"<<G4endl;
132 #endif
133  // parent mass
134  G4double parentmass = current_parent_mass.Get();
135 
136  //create parent G4DynamicParticle at rest
137  G4ThreeVector dummy;
138  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0, parentmass);
139  //create G4Decayproducts
140  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
141  delete parentparticle;
142 
143  //create daughter G4DynamicParticle at rest
144  G4DynamicParticle * daughterparticle = new G4DynamicParticle( G4MT_daughters[0], dummy, 0.0);
145  if (useGivenDaughterMass) daughterparticle->SetMass(givenDaughterMasses[0]);
146  products->PushProducts(daughterparticle);
147 
148 #ifdef G4VERBOSE
149  if (GetVerboseLevel()>1) {
150  G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt ";
151  G4cout << " create decay products in rest frame " <<G4endl;
152  products->DumpInfo();
153  }
154 #endif
155  return products;
156 }
157 
158 G4DecayProducts *G4PhaseSpaceDecayChannel::TwoBodyDecayIt()
159 {
160 #ifdef G4VERBOSE
161  if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt()"<<G4endl;
162 #endif
163  // parent mass
164  G4double parentmass = current_parent_mass.Get();
165 
166  //daughters'mass, width
167  G4double daughtermass[2], daughterwidth[2];
168  daughtermass[0] = G4MT_daughters_mass[0];
169  daughtermass[1] = G4MT_daughters_mass[1];
170  daughterwidth[0] = G4MT_daughters_width[0];
171  daughterwidth[1] = G4MT_daughters_width[1];
172 
173  //create parent G4DynamicParticle at rest
174  G4ThreeVector dummy;
175  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0, parentmass);
176  //create G4Decayproducts
177  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
178  delete parentparticle;
179 
180  if (!useGivenDaughterMass) {
181  G4bool withWidth = (daughterwidth[0]>1.0e-3*daughtermass[0])
182  || (daughterwidth[1]>1.0e-3*daughtermass[1]);
183  if (withWidth) {
184  G4double sumofdaughterwidthsq = daughterwidth[0]*daughterwidth[0]+daughterwidth[1]*daughterwidth[1];
185  // check parent mass and daughter mass
186  G4double maxDev = (parentmass - daughtermass[0] - daughtermass[1] )/std::sqrt(sumofdaughterwidthsq);
187  if (maxDev <= -1.0*rangeMass ){
188 #ifdef G4VERBOSE
189  if (GetVerboseLevel()>0) {
190  G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt "
191  << "sum of daughter mass is larger than parent mass" << G4endl;
192  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass.Get()/GeV << G4endl;
193  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
194  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
195  }
196 #endif
197  G4Exception("G4PhaseSpaceDecayChannel::TwoBodyDecayIt",
198  "PART112", JustWarning,
199  "Can not create decay products: sum of daughter mass is larger than parent mass");
200  return products;
201  }
202  G4double dm1=daughtermass[0];
203  if (daughterwidth[0] > 0.) dm1= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
204  G4double dm2= daughtermass[1];
205  if (daughterwidth[1] > 0.) dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
206  while (dm1+dm2>parentmass){ // Loop checking, 09.08.2015, K.Kurashige
207  dm1= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
208  dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
209  }
210  daughtermass[0] = dm1;
211  daughtermass[1] = dm2;
212  }
213  } else {
214  // use given daughter mass;
215  daughtermass[0] = givenDaughterMasses[0];
216  daughtermass[1] = givenDaughterMasses[1];
217  }
218  if (parentmass < daughtermass[0] + daughtermass[1] ){
219 #ifdef G4VERBOSE
220  if (GetVerboseLevel()>0) {
221  G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt "
222  << "sum of daughter mass is larger than parent mass" << G4endl;
223  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass.Get()/GeV << G4endl;
224  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
225  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
226  if (useGivenDaughterMass) {
227  G4cout << "Daughter Mass is given " << G4endl;
228  }
229  }
230 #endif
231  G4Exception("G4PhaseSpaceDecayChannel::TwoBodyDecayIt",
232  "PART112", JustWarning,
233  "Can not create decay products: sum of daughter mass is larger than parent mass");
234  return products;
235  }
236 
237  //calculate daughter momentum
238  G4double daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
239 
240  G4double costheta = 2.*G4UniformRand()-1.0;
241  G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
242  G4double phi = twopi*G4UniformRand()*rad;
243  G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
244 
245  //create daughter G4DynamicParticle
246  G4double Ekin = std::sqrt(daughtermomentum*daughtermomentum + daughtermass[0]*daughtermass[0]) - daughtermass[0];
247  G4DynamicParticle * daughterparticle = new G4DynamicParticle( G4MT_daughters[0], direction, Ekin, daughtermass[0]);
248  products->PushProducts(daughterparticle);
249  Ekin = std::sqrt(daughtermomentum*daughtermomentum + daughtermass[1]*daughtermass[1]) - daughtermass[1];
250  daughterparticle = new G4DynamicParticle( G4MT_daughters[1], -1.0*direction, Ekin, daughtermass[1]);
251  products->PushProducts(daughterparticle);
252 
253 #ifdef G4VERBOSE
254  if (GetVerboseLevel()>1) {
255  G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt ";
256  G4cout << " create decay products in rest frame " <<G4endl;
257  products->DumpInfo();
258  }
259 #endif
260  return products;
261 }
262 
263 G4DecayProducts *G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()
264 // algorism of this code is originally written in GDECA3 of GEANT3
265 {
266 #ifdef G4VERBOSE
267  if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()"<<G4endl;
268 #endif
269  // parent mass
270  G4double parentmass = current_parent_mass.Get();
271  //daughters'mass
272  G4double daughtermass[3], daughterwidth[3];
273  G4double sumofdaughtermass = 0.0;
274  G4double sumofdaughterwidthsq = 0.0;
275  G4bool withWidth = false;
276  for (G4int index=0; index<3; index++){
277  daughtermass[index] = G4MT_daughters_mass[index];
278  sumofdaughtermass += daughtermass[index];
279  daughterwidth[index] = G4MT_daughters_width[index];
280  sumofdaughterwidthsq += daughterwidth[index]*daughterwidth[index];
281  withWidth = withWidth ||(daughterwidth[index]>1.0e-3*daughtermass[index]);
282  }
283 
284  //create parent G4DynamicParticle at rest
285  G4ThreeVector dummy;
286  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0, parentmass);
287 
288 
289  //create G4Decayproducts
290  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
291  delete parentparticle;
292 
293  if (!useGivenDaughterMass) {
294  if (withWidth){
295  G4double maxDev = (parentmass - sumofdaughtermass )/std::sqrt(sumofdaughterwidthsq) ;
296  if (maxDev <= -1.0*rangeMass ){
297 #ifdef G4VERBOSE
298  if (GetVerboseLevel()>0) {
299  G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt "
300  << "sum of daughter mass is larger than parent mass" << G4endl;
301  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass.Get()/GeV << G4endl;
302  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
303  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
304  G4cout << "daughter 3:" << G4MT_daughters[2]->GetParticleName() << " " << daughtermass[2]/GeV << G4endl;
305  }
306 #endif
307  G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
308  "PART112", JustWarning,
309  "Can not create decay products: sum of daughter mass is larger than parent mass");
310  return products;
311  }
312  G4double dm1=daughtermass[0];
313  if (daughterwidth[0] > 0.) dm1= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
314  G4double dm2= daughtermass[1];
315  if (daughterwidth[1] > 0.) dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
316  G4double dm3= daughtermass[2];
317  if (daughterwidth[2] > 0.) dm3= DynamicalMass(daughtermass[2],daughterwidth[2], maxDev);
318  while (dm1+dm2+dm3>parentmass){ // Loop checking, 09.08.2015, K.Kurashige
319  dm1= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
320  dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
321  dm3= DynamicalMass(daughtermass[2],daughterwidth[2], maxDev);
322  }
323  daughtermass[0] = dm1;
324  daughtermass[1] = dm2;
325  daughtermass[2] = dm3;
326  sumofdaughtermass = dm1+dm2+dm3;
327  }
328  } else {
329  // use given daughter mass;
330  daughtermass[0] = givenDaughterMasses[0];
331  daughtermass[1] = givenDaughterMasses[1];
332  daughtermass[2] = givenDaughterMasses[2];
333  sumofdaughtermass = daughtermass[0] + daughtermass[1] + daughtermass[2];
334  }
335 
336  if (sumofdaughtermass >parentmass) {
337 #ifdef G4VERBOSE
338  if (GetVerboseLevel()>0) {
339  G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt "
340  << "sum of daughter mass is larger than parent mass" << G4endl;
341  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass.Get()/GeV << G4endl;
342  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
343  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
344  G4cout << "daughter 3:" << G4MT_daughters[2]->GetParticleName() << " " << daughtermass[2]/GeV << G4endl;
345  }
346  if (useGivenDaughterMass) {
347  G4cout << "Daughter Mass is given " << G4endl;
348  }
349 #endif
350  G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
351  "PART112", JustWarning,
352  "Can not create decay products: sum of daughter mass is larger than parent mass");
353  return products;
354  }
355 
356 
357 
358  //calculate daughter momentum
359  // Generate two
360  G4double rd1, rd2, rd;
361  G4double daughtermomentum[3];
362  G4double momentummax=0.0, momentumsum = 0.0;
363  G4double energy;
364  const size_t MAX_LOOP=10000;
365 
366  for (size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
367  rd1 = G4UniformRand();
368  rd2 = G4UniformRand();
369  if (rd2 > rd1) {
370  rd = rd1;
371  rd1 = rd2;
372  rd2 = rd;
373  }
374  momentummax = 0.0;
375  momentumsum = 0.0;
376  // daughter 0
377  energy = rd2*(parentmass - sumofdaughtermass);
378  daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
379  if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0];
380  momentumsum += daughtermomentum[0];
381  // daughter 1
382  energy = (1.-rd1)*(parentmass - sumofdaughtermass);
383  daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
384  if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1];
385  momentumsum += daughtermomentum[1];
386  // daughter 2
387  energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
388  daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
389  if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2];
390  momentumsum += daughtermomentum[2];
391  if (momentummax <= momentumsum - momentummax ) break;
392  }
393 
394  // output message
395 #ifdef G4VERBOSE
396  if (GetVerboseLevel()>1) {
397  G4cout << " daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl;
398  G4cout << " daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl;
399  G4cout << " daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl;
400  G4cout << " momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl;
401  }
402 #endif
403 
404  //create daughter G4DynamicParticle
405  G4double costheta, sintheta, phi, sinphi, cosphi;
406  G4double costhetan, sinthetan, phin, sinphin, cosphin;
407  costheta = 2.*G4UniformRand()-1.0;
408  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
409  phi = twopi*G4UniformRand()*rad;
410  sinphi = std::sin(phi);
411  cosphi = std::cos(phi);
412 
413  G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
414  G4double Ekin = std::sqrt(daughtermomentum[0]*daughtermomentum[0] + daughtermass[0]*daughtermass[0]) - daughtermass[0];
415  G4DynamicParticle * daughterparticle = new G4DynamicParticle( G4MT_daughters[0],
416  direction0,
417  Ekin, daughtermass[0]);
418  products->PushProducts(daughterparticle);
419 
420  costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
421  sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
422  phin = twopi*G4UniformRand()*rad;
423  sinphin = std::sin(phin);
424  cosphin = std::cos(phin);
425  G4ThreeVector direction2;
426  direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
427  direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
428  direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
429  G4ThreeVector pmom = daughtermomentum[2]*direction2/direction2.mag();
430  Ekin = std::sqrt(pmom.mag2() + daughtermass[2]*daughtermass[2]) - daughtermass[2];
431  daughterparticle = new G4DynamicParticle( G4MT_daughters[2],
432  pmom/pmom.mag(),
433  Ekin, daughtermass[2]);
434  products->PushProducts(daughterparticle);
435 
436  pmom = (direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0);
437  Ekin = std::sqrt(pmom.mag2() + daughtermass[1]*daughtermass[1]) - daughtermass[1];
438  daughterparticle = new G4DynamicParticle(
439  G4MT_daughters[1],
440  pmom/pmom.mag(),
441  Ekin, daughtermass[1]);
442  products->PushProducts(daughterparticle);
443 
444 #ifdef G4VERBOSE
445  if (GetVerboseLevel()>1) {
446  G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt ";
447  G4cout << " create decay products in rest frame " <<G4endl;
448  products->DumpInfo();
449  }
450 #endif
451  return products;
452 }
453 
454 G4DecayProducts *G4PhaseSpaceDecayChannel::ManyBodyDecayIt()
455 // algorithm of this code is originally written in FORTRAN by M.Asai
456 //******************************************************************
457 // NBODY
458 // N-body phase space Monte-Carlo generator
459 // Makoto Asai
460 // Hiroshima Institute of Technology
461 // (asai@kekvax.kek.jp)
462 // Revised release : 19/Apr/1995
463 //
464 {
465  G4int index, index2;
466 
467 #ifdef G4VERBOSE
468  if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()"<<G4endl;
469 #endif
470 
471 
472  // parent mass
473  G4double parentmass = current_parent_mass.Get();
474 
475  //parent particle
476  G4ThreeVector dummy;
477  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0, parentmass);
478 
479  //create G4Decayproducts
480  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
481  delete parentparticle;
482 
483  //daughters'mass
484  G4double *daughtermass = new G4double[numberOfDaughters];
485 
486  G4double sumofdaughtermass = 0.0;
487  for (index=0; index<numberOfDaughters; index++){
488  if (!useGivenDaughterMass) {
489  daughtermass[index] = G4MT_daughters_mass[index];
490  } else {
491  daughtermass[index] = givenDaughterMasses[index];
492  }
493  sumofdaughtermass += daughtermass[index];
494  }
495  if (sumofdaughtermass >parentmass) {
496 #ifdef G4VERBOSE
497  if (GetVerboseLevel()>0) {
498  G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt "
499  << "sum of daughter mass is larger than parent mass" << G4endl;
500  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass.Get()/GeV << G4endl;
501  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
502  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
503  G4cout << "daughter 3:" << G4MT_daughters[2]->GetParticleName() << " " << daughtermass[2]/GeV << G4endl;
504  G4cout << "daughter 4:" << G4MT_daughters[3]->GetParticleName() << " " << daughtermass[3]/GeV << G4endl;
505  }
506 #endif
507  G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt",
508  "PART112", JustWarning,
509  "Can not create decay products: sum of daughter mass is larger than parent mass");
510  delete [] daughtermass;
511  return products;
512  }
513 
514  //Calculate daughter momentum
515  G4double *daughtermomentum = new G4double[numberOfDaughters];
516  G4ThreeVector direction;
517  G4DynamicParticle **daughterparticle;
518  G4double *sm = new G4double[numberOfDaughters];
519  G4double tmas;
520  G4double weight = 1.0;
521  G4int numberOfTry = 0;
522  const size_t MAX_LOOP=10000;
523 
524  for (size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
525  //Generate rundom number in descending order
526  G4double temp;
527  G4double *rd = new G4double[numberOfDaughters];
528  rd[0] = 1.0;
529  for(index =1; index < numberOfDaughters -1; index++) rd[index] = G4UniformRand();
530  rd[ numberOfDaughters -1] = 0.0;
531  for(index =1; index < numberOfDaughters -1; index++) {
532  for(index2 = index+1; index2 < numberOfDaughters; index2++) {
533  if (rd[index] < rd[index2]){
534  temp = rd[index];
535  rd[index] = rd[index2];
536  rd[index2] = temp;
537  }
538  }
539  }
540 
541  //calcurate virtual mass
542  tmas = parentmass - sumofdaughtermass;
543  temp = sumofdaughtermass;
544  for(index =0; index < numberOfDaughters; index++) {
545  sm[index] = rd[index]*tmas + temp;
546  temp -= daughtermass[index];
547  if (GetVerboseLevel()>1) {
548  G4cout << index << " rundom number:" << rd[index];
549  G4cout << " virtual mass:" << sm[index]/GeV << "[GeV/c/c]" <<G4endl;
550  }
551  }
552  delete [] rd;
553 
554  //Calculate daughter momentum
555  weight = 1.0;
556  G4bool smOK=true;
557  for(index =0; index< numberOfDaughters-1 && smOK; index--) {
558  smOK = (sm[index]-daughtermass[index]- sm[index+1] >=0.);
559  }
560  if (!smOK) continue;
561 
562  index =numberOfDaughters-1;
563  daughtermomentum[index]= Pmx( sm[index-1],daughtermass[index-1], sm[index]);
564 #ifdef G4VERBOSE
565  if (GetVerboseLevel()>1) {
566  G4cout << " daughter " << index << ":" << *daughters_name[index];
567  G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
568  }
569 #endif
570  for(index =numberOfDaughters-2; index>=0; index--) {
571  // calculate
572  daughtermomentum[index]= Pmx( sm[index],daughtermass[index], sm[index +1]);
573  if(daughtermomentum[index] < 0.0) {
574  // !!! illegal momentum !!!
575 #ifdef G4VERBOSE
576  if (GetVerboseLevel()>0) {
577  G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
578  G4cout << " can not calculate daughter momentum " <<G4endl;
579  G4cout << " parent:" << *parent_name;
580  G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
581  G4cout << " daughter " << index << ":" << *daughters_name[index];
582  G4cout << " mass:" << daughtermass[index]/GeV << "[GeV/c/c]" ;
583  G4cout << " mass:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
584  if (useGivenDaughterMass) {
585  G4cout << "Daughter Mass is given " << G4endl;
586  }
587  }
588 #endif
589  delete [] sm;
590  delete [] daughtermass;
591  delete [] daughtermomentum;
592  delete products;
593  G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt",
594  "PART112", JustWarning,
595  "Can not create decay products: sum of daughter mass is larger than parent mass");
596 
597  return 0; // Error detection
598 
599  } else {
600  // calculate weight of this events
601  weight *= daughtermomentum[index]/sm[index];
602 #ifdef G4VERBOSE
603  if (GetVerboseLevel()>1) {
604  G4cout << " daughter " << index << ":" << *daughters_name[index];
605  G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
606  }
607 #endif
608  }
609  }
610 
611 #ifdef G4VERBOSE
612  if (GetVerboseLevel()>1) {
613  G4cout << " weight: " << weight <<G4endl;
614  }
615 #endif
616 
617  // exit if number of Try exceeds 100
618  if (numberOfTry++ >100) {
619 #ifdef G4VERBOSE
620  if (GetVerboseLevel()>0) {
621  G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ";
622  G4cout << " can not determine Decay Kinematics " << G4endl;
623  G4cout << "parent : " << *parent_name << G4endl;
624  G4cout << "daughters : ";
625  for(index=0; index<numberOfDaughters; index++) {
626  G4cout << *daughters_name[index] << " , ";
627  }
628  G4cout << G4endl;
629  }
630 #endif
631  G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ",
632  "PART113", JustWarning,
633  " Cannot decay : Decay Kinematics cannot be calculated ");
634 
635  delete [] sm;
636  delete [] daughtermass;
637  delete [] daughtermomentum;
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 
G4PhaseSpaceDecayChannel::DecayIt
virtual G4DecayProducts * DecayIt(G4double)
Definition: G4PhaseSpaceDecayChannel.cc:82
G4LorentzVector.hh
G4VDecayChannel::CheckAndFillDaughters
void CheckAndFillDaughters()
Definition: G4VDecayChannel.hh:324
G4VDecayChannel::G4MT_daughters_mass
G4double * G4MT_daughters_mass
Definition: G4VDecayChannel.hh:195
G4PhaseSpaceDecayChannel::SetDaughterMasses
G4bool SetDaughterMasses(G4double masses[])
Definition: G4PhaseSpaceDecayChannel.cc:711
CLHEP::Hep3Vector::x
double x() const
G4DecayProducts::PushProducts
G4int PushProducts(G4DynamicParticle *aParticle)
Definition: G4DecayProducts.cc:153
G4VDecayChannel::G4MT_parent
G4ParticleDefinition * G4MT_parent
Definition: G4VDecayChannel.hh:192
G4PhaseSpaceDecayChannel::IsOKWithParentMass
G4bool IsOKWithParentMass(G4double parentMass)
Definition: G4PhaseSpaceDecayChannel.cc:726
G4Cache::Get
value_type & Get() const
Definition: G4Cache.hh:282
G4VDecayChannel::G4MT_daughters
G4ParticleDefinition ** G4MT_daughters
Definition: G4VDecayChannel.hh:193
G4VDecayChannel::G4MT_parent_mass
G4double G4MT_parent_mass
Definition: G4VDecayChannel.hh:194
G4VDecayChannel::IsOKWithParentMass
virtual G4bool IsOKWithParentMass(G4double parentMass)
Definition: G4VDecayChannel.cc:586
rad
static constexpr double rad
Definition: G4SIunits.hh:149
JustWarning
Definition: G4ExceptionSeverity.hh:64
G4int
int G4int
Definition: G4Types.hh:78
CLHEP::Hep3Vector::setY
void setY(double)
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:120
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::setZ
void setZ(double)
CLHEP::Hep3Vector::setX
void setX(double)
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
G4VDecayChannel.hh
G4bool
bool G4bool
Definition: G4Types.hh:79
G4DecayProducts::DumpInfo
void DumpInfo() const
Definition: G4DecayProducts.cc:311
CLHEP::HepLorentzVector::boost
HepLorentzVector & boost(double, double, double)
Definition: LorentzVector.cc:59
G4VDecayChannel::parent_name
G4String * parent_name
Definition: G4VDecayChannel.hh:171
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
G4DecayProducts.hh
G4DynamicParticle::Set4Momentum
void Set4Momentum(const G4LorentzVector &momentum)
Definition: G4DynamicParticle.cc:385
G4VDecayChannel::GetVerboseLevel
G4int GetVerboseLevel() const
Definition: G4VDecayChannel.hh:303
G4VDecayChannel::G4MT_daughters_width
G4double * G4MT_daughters_width
Definition: G4VDecayChannel.hh:196
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:157
CLHEP::Hep3Vector::y
double y() const
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
CLHEP::Hep3Vector::mag2
double mag2() const
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VDecayChannel::DynamicalMass
G4double DynamicalMass(G4double massPDG, G4double width, G4double maxDev=+1.) const
Definition: G4VDecayChannel.cc:568
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
CLHEP::Hep3Vector::mag
double mag() const
G4Cache::Put
void Put(const value_type &val) const
Definition: G4Cache.hh:286
G4DynamicParticle::SetMass
void SetMass(G4double mass)
Definition: G4DynamicParticle.cc:445
G4VDecayChannel::daughters_name
G4String ** daughters_name
Definition: G4VDecayChannel.hh:173
G4PhaseSpaceDecayChannel::Pmx
static G4double Pmx(G4double e, G4double p1, G4double p2)
Definition: G4PhaseSpaceDecayChannel.cc:739
G4LorentzRotation.hh