Geant4  10.02.p02
G4PhaseSpaceDecayChannel.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4PhaseSpaceDecayChannel.cc 97537 2016-06-03 15:26:56Z 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  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= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
202  G4double dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
203  while (dm1+dm2>parentmass){ // Loop checking, 09.08.2015, K.Kurashige
204  dm1= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
205  dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
206  }
207  daughtermass[0] = dm1;
208  daughtermass[1] = dm2;
209  }
210  } else {
211  // use given daughter mass;
212  daughtermass[0] = givenDaughterMasses[0];
213  daughtermass[1] = givenDaughterMasses[1];
214  }
215  if (parentmass < daughtermass[0] + daughtermass[1] ){
216 #ifdef G4VERBOSE
217  if (GetVerboseLevel()>0) {
218  G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt "
219  << "sum of daughter mass is larger than parent mass" << G4endl;
220  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass.Get()/GeV << G4endl;
221  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
222  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
223  if (useGivenDaughterMass) {
224  G4cout << "Daughter Mass is given " << G4endl;
225  }
226  }
227 #endif
228  G4Exception("G4PhaseSpaceDecayChannel::TwoBodyDecayIt",
229  "PART112", JustWarning,
230  "Can not create decay products: sum of daughter mass is larger than parent mass");
231  return products;
232  }
233 
234  //calculate daughter momentum
235  G4double daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
236 
237  G4double costheta = 2.*G4UniformRand()-1.0;
238  G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
239  G4double phi = twopi*G4UniformRand()*rad;
240  G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
241 
242  //create daughter G4DynamicParticle
243  G4double Ekin = std::sqrt(daughtermomentum*daughtermomentum + daughtermass[0]*daughtermass[0]) - daughtermass[0];
244  G4DynamicParticle * daughterparticle = new G4DynamicParticle( G4MT_daughters[0], direction, Ekin, daughtermass[0]);
245  products->PushProducts(daughterparticle);
246  Ekin = std::sqrt(daughtermomentum*daughtermomentum + daughtermass[1]*daughtermass[1]) - daughtermass[1];
247  daughterparticle = new G4DynamicParticle( G4MT_daughters[1], -1.0*direction, Ekin, daughtermass[1]);
248  products->PushProducts(daughterparticle);
249 
250 #ifdef G4VERBOSE
251  if (GetVerboseLevel()>1) {
252  G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt ";
253  G4cout << " create decay products in rest frame " <<G4endl;
254  products->DumpInfo();
255  }
256 #endif
257  return products;
258 }
259 
260 G4DecayProducts *G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()
261 // algorism of this code is originally written in GDECA3 of GEANT3
262 {
263 #ifdef G4VERBOSE
264  if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()"<<G4endl;
265 #endif
266  // parent mass
267  G4double parentmass = current_parent_mass.Get();
268  //daughters'mass
269  G4double daughtermass[3], daughterwidth[3];
270  G4double sumofdaughtermass = 0.0;
271  G4double sumofdaughterwidthsq = 0.0;
272  G4bool withWidth = false;
273  for (G4int index=0; index<3; index++){
274  daughtermass[index] = G4MT_daughters_mass[index];
275  sumofdaughtermass += daughtermass[index];
276  daughterwidth[index] = G4MT_daughters_width[index];
277  sumofdaughterwidthsq += daughterwidth[index]*daughterwidth[index];
278  withWidth = withWidth ||(daughterwidth[index]>1.0e-3*daughtermass[index]);
279  }
280 
281  //create parent G4DynamicParticle at rest
282  G4ThreeVector dummy;
283  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0, parentmass);
284 
285 
286  //create G4Decayproducts
287  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
288  delete parentparticle;
289 
290  if (!useGivenDaughterMass) {
291  if (withWidth){
292  G4double maxDev = (parentmass - sumofdaughtermass )/std::sqrt(sumofdaughterwidthsq) ;
293  if (maxDev <= -1.0*rangeMass ){
294 #ifdef G4VERBOSE
295  if (GetVerboseLevel()>0) {
296  G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt "
297  << "sum of daughter mass is larger than parent mass" << G4endl;
298  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass.Get()/GeV << G4endl;
299  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
300  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
301  G4cout << "daughter 3:" << G4MT_daughters[2]->GetParticleName() << " " << daughtermass[2]/GeV << G4endl;
302  }
303 #endif
304  G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
305  "PART112", JustWarning,
306  "Can not create decay products: sum of daughter mass is larger than parent mass");
307  return products;
308  }
309  G4double dm1= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
310  G4double dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
311  G4double dm3= DynamicalMass(daughtermass[2],daughterwidth[2], maxDev);
312  while (dm1+dm2+dm3>parentmass){ // Loop checking, 09.08.2015, K.Kurashige
313  dm1= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
314  dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
315  dm3= DynamicalMass(daughtermass[2],daughterwidth[2], maxDev);
316  }
317  daughtermass[0] = dm1;
318  daughtermass[1] = dm2;
319  daughtermass[2] = dm3;
320  sumofdaughtermass = dm1+dm2+dm3;
321  }
322  } else {
323  // use given daughter mass;
324  daughtermass[0] = givenDaughterMasses[0];
325  daughtermass[1] = givenDaughterMasses[1];
326  daughtermass[2] = givenDaughterMasses[2];
327  sumofdaughtermass = daughtermass[0] + daughtermass[1] + daughtermass[2];
328  }
329 
330  if (sumofdaughtermass >parentmass) {
331 #ifdef G4VERBOSE
332  if (GetVerboseLevel()>0) {
333  G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt "
334  << "sum of daughter mass is larger than parent mass" << G4endl;
335  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass.Get()/GeV << G4endl;
336  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
337  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
338  G4cout << "daughter 3:" << G4MT_daughters[2]->GetParticleName() << " " << daughtermass[2]/GeV << G4endl;
339  }
340  if (useGivenDaughterMass) {
341  G4cout << "Daughter Mass is given " << G4endl;
342  }
343 #endif
344  G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
345  "PART112", JustWarning,
346  "Can not create decay products: sum of daughter mass is larger than parent mass");
347  return products;
348  }
349 
350 
351 
352  //calculate daughter momentum
353  // Generate two
354  G4double rd1, rd2, rd;
355  G4double daughtermomentum[3];
356  G4double momentummax=0.0, momentumsum = 0.0;
357  G4double energy;
358  const size_t MAX_LOOP=10000;
359 
360  for (size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
361  rd1 = G4UniformRand();
362  rd2 = G4UniformRand();
363  if (rd2 > rd1) {
364  rd = rd1;
365  rd1 = rd2;
366  rd2 = rd;
367  }
368  momentummax = 0.0;
369  momentumsum = 0.0;
370  // daughter 0
371  energy = rd2*(parentmass - sumofdaughtermass);
372  daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
373  if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0];
374  momentumsum += daughtermomentum[0];
375  // daughter 1
376  energy = (1.-rd1)*(parentmass - sumofdaughtermass);
377  daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
378  if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1];
379  momentumsum += daughtermomentum[1];
380  // daughter 2
381  energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
382  daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
383  if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2];
384  momentumsum += daughtermomentum[2];
385  if (momentummax <= momentumsum - momentummax ) break;
386  }
387 
388  // output message
389 #ifdef G4VERBOSE
390  if (GetVerboseLevel()>1) {
391  G4cout << " daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl;
392  G4cout << " daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl;
393  G4cout << " daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl;
394  G4cout << " momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl;
395  }
396 #endif
397 
398  //create daughter G4DynamicParticle
399  G4double costheta, sintheta, phi, sinphi, cosphi;
400  G4double costhetan, sinthetan, phin, sinphin, cosphin;
401  costheta = 2.*G4UniformRand()-1.0;
402  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
403  phi = twopi*G4UniformRand()*rad;
404  sinphi = std::sin(phi);
405  cosphi = std::cos(phi);
406 
407  G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
408  G4double Ekin = std::sqrt(daughtermomentum[0]*daughtermomentum[0] + daughtermass[0]*daughtermass[0]) - daughtermass[0];
409  G4DynamicParticle * daughterparticle = new G4DynamicParticle( G4MT_daughters[0],
410  direction0,
411  Ekin, daughtermass[0]);
412  products->PushProducts(daughterparticle);
413 
414  costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
415  sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
416  phin = twopi*G4UniformRand()*rad;
417  sinphin = std::sin(phin);
418  cosphin = std::cos(phin);
419  G4ThreeVector direction2;
420  direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
421  direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
422  direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
423  G4ThreeVector pmom = daughtermomentum[2]*direction2/direction2.mag();
424  Ekin = std::sqrt(pmom.mag2() + daughtermass[2]*daughtermass[2]) - daughtermass[2];
425  daughterparticle = new G4DynamicParticle( G4MT_daughters[2],
426  pmom/pmom.mag(),
427  Ekin, daughtermass[2]);
428  products->PushProducts(daughterparticle);
429 
430  pmom = (direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0);
431  Ekin = std::sqrt(pmom.mag2() + daughtermass[1]*daughtermass[1]) - daughtermass[1];
432  daughterparticle = new G4DynamicParticle(
433  G4MT_daughters[1],
434  pmom/pmom.mag(),
435  Ekin, daughtermass[1]);
436  products->PushProducts(daughterparticle);
437 
438 #ifdef G4VERBOSE
439  if (GetVerboseLevel()>1) {
440  G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt ";
441  G4cout << " create decay products in rest frame " <<G4endl;
442  products->DumpInfo();
443  }
444 #endif
445  return products;
446 }
447 
448 G4DecayProducts *G4PhaseSpaceDecayChannel::ManyBodyDecayIt()
449 // algorithm of this code is originally written in FORTRAN by M.Asai
450 //******************************************************************
451 // NBODY
452 // N-body phase space Monte-Carlo generator
453 // Makoto Asai
454 // Hiroshima Institute of Technology
455 // (asai@kekvax.kek.jp)
456 // Revised release : 19/Apr/1995
457 //
458 {
459  G4int index, index2;
460 
461 #ifdef G4VERBOSE
462  if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()"<<G4endl;
463 #endif
464 
465 
466  // parent mass
467  G4double parentmass = current_parent_mass.Get();
468 
469  //parent particle
470  G4ThreeVector dummy;
471  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0, parentmass);
472 
473  //create G4Decayproducts
474  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
475  delete parentparticle;
476 
477  //daughters'mass
478  G4double *daughtermass = new G4double[numberOfDaughters];
479 
480  G4double sumofdaughtermass = 0.0;
481  for (index=0; index<numberOfDaughters; index++){
482  if (!useGivenDaughterMass) {
483  daughtermass[index] = G4MT_daughters_mass[index];
484  } else {
485  daughtermass[index] = givenDaughterMasses[index];
486  }
487  sumofdaughtermass += daughtermass[index];
488  }
489  if (sumofdaughtermass >parentmass) {
490 #ifdef G4VERBOSE
491  if (GetVerboseLevel()>0) {
492  G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt "
493  << "sum of daughter mass is larger than parent mass" << G4endl;
494  G4cout << "parent :" << G4MT_parent->GetParticleName() << " " << current_parent_mass.Get()/GeV << G4endl;
495  G4cout << "daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
496  G4cout << "daughter 2:" << G4MT_daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
497  G4cout << "daughter 3:" << G4MT_daughters[2]->GetParticleName() << " " << daughtermass[2]/GeV << G4endl;
498  G4cout << "daughter 4:" << G4MT_daughters[3]->GetParticleName() << " " << daughtermass[3]/GeV << G4endl;
499  }
500 #endif
501  G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt",
502  "PART112", JustWarning,
503  "Can not create decay products: sum of daughter mass is larger than parent mass");
504  return products;
505  }
506 
507  //Calculate daughter momentum
508  G4double *daughtermomentum = new G4double[numberOfDaughters];
509  G4ThreeVector direction;
510  G4DynamicParticle **daughterparticle;
511  G4double *sm = new G4double[numberOfDaughters];
512  G4double tmas;
513  G4double weight = 1.0;
514  G4int numberOfTry = 0;
515  const size_t MAX_LOOP=10000;
516 
517  for (size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
518  //Generate rundom number in descending order
519  G4double temp;
520  G4double *rd = new G4double[numberOfDaughters];
521  rd[0] = 1.0;
522  for(index =1; index < numberOfDaughters -1; index++) rd[index] = G4UniformRand();
523  rd[ numberOfDaughters -1] = 0.0;
524  for(index =1; index < numberOfDaughters -1; index++) {
525  for(index2 = index+1; index2 < numberOfDaughters; index2++) {
526  if (rd[index] < rd[index2]){
527  temp = rd[index];
528  rd[index] = rd[index2];
529  rd[index2] = temp;
530  }
531  }
532  }
533 
534  //calcurate virtual mass
535  tmas = parentmass - sumofdaughtermass;
536  temp = sumofdaughtermass;
537  for(index =0; index < numberOfDaughters; index++) {
538  sm[index] = rd[index]*tmas + temp;
539  temp -= daughtermass[index];
540  if (GetVerboseLevel()>1) {
541  G4cout << index << " rundom number:" << rd[index];
542  G4cout << " virtual mass:" << sm[index]/GeV << "[GeV/c/c]" <<G4endl;
543  }
544  }
545  delete [] rd;
546 
547  //Calculate daughter momentum
548  weight = 1.0;
549  G4bool smOK=true;
550  for(index =0; index< numberOfDaughters-1 && smOK; index--) {
551  smOK = (sm[index]-daughtermass[index]- sm[index+1] >=0.);
552  }
553  if (!smOK) continue;
554 
555  index =numberOfDaughters-1;
556  daughtermomentum[index]= Pmx( sm[index-1],daughtermass[index-1], sm[index]);
557 #ifdef G4VERBOSE
558  if (GetVerboseLevel()>1) {
559  G4cout << " daughter " << index << ":" << *daughters_name[index];
560  G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
561  }
562 #endif
563  for(index =numberOfDaughters-2; index>=0; index--) {
564  // calculate
565  daughtermomentum[index]= Pmx( sm[index],daughtermass[index], sm[index +1]);
566  if(daughtermomentum[index] < 0.0) {
567  // !!! illegal momentum !!!
568 #ifdef G4VERBOSE
569  if (GetVerboseLevel()>0) {
570  G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
571  G4cout << " can not calculate daughter momentum " <<G4endl;
572  G4cout << " parent:" << *parent_name;
573  G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
574  G4cout << " daughter " << index << ":" << *daughters_name[index];
575  G4cout << " mass:" << daughtermass[index]/GeV << "[GeV/c/c]" ;
576  G4cout << " mass:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
577  if (useGivenDaughterMass) {
578  G4cout << "Daughter Mass is given " << G4endl;
579  }
580  }
581 #endif
582  delete [] sm;
583  delete [] daughtermass;
584  delete [] daughtermomentum;
585 
586  G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt",
587  "PART112", JustWarning,
588  "Can not create decay products: sum of daughter mass is larger than parent mass");
589 
590  return 0; // Error detection
591 
592  } else {
593  // calculate weight of this events
594  weight *= daughtermomentum[index]/sm[index];
595 #ifdef G4VERBOSE
596  if (GetVerboseLevel()>1) {
597  G4cout << " daughter " << index << ":" << *daughters_name[index];
598  G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
599  }
600 #endif
601  }
602  }
603 
604 #ifdef G4VERBOSE
605  if (GetVerboseLevel()>1) {
606  G4cout << " weight: " << weight <<G4endl;
607  }
608 #endif
609 
610  // exit if number of Try exceeds 100
611  if (numberOfTry++ >100) {
612 #ifdef G4VERBOSE
613  if (GetVerboseLevel()>0) {
614  G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ";
615  G4cout << " can not determine Decay Kinematics " << G4endl;
616  G4cout << "parent : " << *parent_name << G4endl;
617  G4cout << "daughters : ";
618  for(index=0; index<numberOfDaughters; index++) {
619  G4cout << *daughters_name[index] << " , ";
620  }
621  G4cout << G4endl;
622  }
623 #endif
624  G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ",
625  "PART113", JustWarning,
626  " Cannot decay : Decay Kinematics cannot be calculated ");
627 
628  delete [] sm;
629  delete [] daughtermass;
630  delete [] daughtermomentum;
631  return 0; // Error detection
632  }
633  if ( weight < G4UniformRand()) break;
634  }
635 
636 #ifdef G4VERBOSE
637  if (GetVerboseLevel()>1) {
638  G4cout << "Start calulation of daughters momentum vector "<<G4endl;
639  }
640 #endif
641 
642  G4double costheta, sintheta, phi;
643  G4double beta;
644  daughterparticle = new G4DynamicParticle*[numberOfDaughters];
645 
646  index = numberOfDaughters -2;
647  costheta = 2.*G4UniformRand()-1.0;
648  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
649  phi = twopi*G4UniformRand()*rad;
650  direction.setZ(costheta);
651  direction.setY(sintheta*std::sin(phi));
652  direction.setX(sintheta*std::cos(phi));
653  daughterparticle[index] = new G4DynamicParticle( G4MT_daughters[index], direction*daughtermomentum[index] );
654  daughterparticle[index+1] = new G4DynamicParticle( G4MT_daughters[index+1], direction*(-1.0*daughtermomentum[index]) );
655 
656  for (index = numberOfDaughters -3; index >= 0; index--) {
657  //calculate momentum direction
658  costheta = 2.*G4UniformRand()-1.0;
659  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
660  phi = twopi*G4UniformRand()*rad;
661  direction.setZ(costheta);
662  direction.setY(sintheta*std::sin(phi));
663  direction.setX(sintheta*std::cos(phi));
664 
665  // boost already created particles
666  beta = daughtermomentum[index];
667  beta /= std::sqrt( daughtermomentum[index]*daughtermomentum[index] + sm[index+1]*sm[index+1] );
668  for (index2 = index+1; index2<numberOfDaughters; index2++) {
669  G4LorentzVector p4;
670  // make G4LorentzVector for secondaries
671  p4 = daughterparticle[index2]->Get4Momentum();
672 
673  // boost secondaries to new frame
674  p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
675 
676  // change energy/momentum
677  daughterparticle[index2]->Set4Momentum(p4);
678  }
679  //create daughter G4DynamicParticle
680  daughterparticle[index]= new G4DynamicParticle( G4MT_daughters[index], direction*(-1.0*daughtermomentum[index]));
681  }
682 
683  //add daughters to G4Decayproducts
684  for (index = 0; index<numberOfDaughters; index++) {
685  products->PushProducts(daughterparticle[index]);
686  }
687 
688 #ifdef G4VERBOSE
689  if (GetVerboseLevel()>1) {
690  G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
691  G4cout << " create decay products in rest frame " << G4endl;
692  products->DumpInfo();
693  }
694 #endif
695 
696  delete [] daughterparticle;
697  delete [] daughtermomentum;
698  delete [] daughtermass;
699  delete [] sm;
700 
701  return products;
702 }
703 
704 G4bool G4PhaseSpaceDecayChannel::SetDaughterMasses( G4double masses[])
705 {
706  for (G4int idx=0; idx<numberOfDaughters; idx++){
707  givenDaughterMasses[idx] = masses[idx];
708  }
709  useGivenDaughterMass = true;
710  return useGivenDaughterMass;
711 }
712 
713 G4bool G4PhaseSpaceDecayChannel::SampleDaughterMasses()
714 {
715  useGivenDaughterMass = false;
716  return useGivenDaughterMass;
717 }
718 
719 G4bool G4PhaseSpaceDecayChannel::IsOKWithParentMass(G4double parentMass) {
720  if (!useGivenDaughterMass) return G4VDecayChannel::IsOKWithParentMass(parentMass);
721 
722  CheckAndFillParent();
723  CheckAndFillDaughters();
724 
725  G4double sumOfDaughterMassMin=0.0;
726  for (G4int index=0; index < numberOfDaughters; index++) {
727  sumOfDaughterMassMin += givenDaughterMasses[index];
728  }
729  return (parentMass >= sumOfDaughterMassMin);
730 }
731 
732 G4double G4PhaseSpaceDecayChannel::Pmx(G4double e, G4double p1, G4double p2)
733 {
734  // calcurate momentum of daughter particles in two-body decay
735  G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e);
736  if (ppp>0) return std::sqrt(ppp);
737  else return -1.;
738 }
739 
740 
741 
G4PhaseSpaceDecayChannel::DecayIt
virtual G4DecayProducts * DecayIt(G4double)
Definition: G4PhaseSpaceDecayChannel.cc:81
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:704
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:85
G4VDecayChannel::G4MT_parent
G4ParticleDefinition * G4MT_parent
Definition: G4VDecayChannel.hh:192
G4PhaseSpaceDecayChannel::current_parent_mass
G4Cache< G4double > current_parent_mass
Definition: G4PhaseSpaceDecayChannel.hh:79
G4PhaseSpaceDecayChannel::IsOKWithParentMass
G4bool IsOKWithParentMass(G4double parentMass)
Definition: G4PhaseSpaceDecayChannel.cc:719
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:584
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:171
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
G4DynamicParticle::Set4Momentum
void Set4Momentum(const G4LorentzVector &momentum)
Definition: G4DynamicParticle.cc:376
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
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VDecayChannel::DynamicalMass
G4double DynamicalMass(G4double massPDG, G4double width, G4double maxDev=+1.) const
Definition: G4VDecayChannel.cc:567
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4Cache::Put
void Put(const value_type &val) const
Definition: G4Cache.hh:286
G4DynamicParticle::SetMass
void SetMass(G4double mass)
G4VDecayChannel::daughters_name
G4String ** daughters_name
Definition: G4VDecayChannel.hh:173
G4PhaseSpaceDecayChannel::Pmx
static G4double Pmx(G4double e, G4double p1, G4double p2)
Definition: G4PhaseSpaceDecayChannel.cc:732
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
G4LorentzRotation.hh