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