ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GeneralPhaseSpaceDecay.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4GeneralPhaseSpaceDecay.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // ----------------------------------------------------------------
28 // GEANT 4 class header file
29 //
30 // History: first implementation, A. Feliciello, 21st May 1998
31 //
32 // Note: this class is a generalization of the
33 // G4PhaseSpaceDecayChannel one
34 // ----------------------------------------------------------------
35 
36 #include "G4ParticleDefinition.hh"
37 #include "G4DecayProducts.hh"
38 #include "G4VDecayChannel.hh"
39 #include "G4GeneralPhaseSpaceDecay.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "Randomize.hh"
43 #include "G4LorentzVector.hh"
44 #include "G4LorentzRotation.hh"
45 #include "G4ios.hh"
46 
47 
48 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay(G4int Verbose) :
49  G4VDecayChannel("Phase Space", Verbose),
50  parentmass(0.), theDaughterMasses(0)
51 {
52  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
53 }
54 
55 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay(const G4String& theParentName,
56  G4double theBR,
57  G4int theNumberOfDaughters,
58  const G4String& theDaughterName1,
59  const G4String& theDaughterName2,
60  const G4String& theDaughterName3) :
61  G4VDecayChannel("Phase Space",
62  theParentName,theBR,
63  theNumberOfDaughters,
64  theDaughterName1,
65  theDaughterName2,
66  theDaughterName3),
67  theDaughterMasses(0)
68 {
69  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
70 
71  // Set the parent particle (resonance) mass to the (default) PDG vale
72  if (G4MT_parent != NULL)
73  {
74  parentmass = G4MT_parent->GetPDGMass();
75  } else {
76  parentmass=0.;
77  }
78 
79 }
80 
81 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay(const G4String& theParentName,
82  G4double theParentMass,
83  G4double theBR,
84  G4int theNumberOfDaughters,
85  const G4String& theDaughterName1,
86  const G4String& theDaughterName2,
87  const G4String& theDaughterName3) :
88  G4VDecayChannel("Phase Space",
89  theParentName,theBR,
90  theNumberOfDaughters,
91  theDaughterName1,
92  theDaughterName2,
93  theDaughterName3),
94  parentmass(theParentMass),
95  theDaughterMasses(0)
96 {
97  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
98 }
99 
100 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay(const G4String& theParentName,
101  G4double theParentMass,
102  G4double theBR,
103  G4int theNumberOfDaughters,
104  const G4String& theDaughterName1,
105  const G4String& theDaughterName2,
106  const G4String& theDaughterName3,
107  const G4double *masses) :
108  G4VDecayChannel("Phase Space",
109  theParentName,theBR,
110  theNumberOfDaughters,
111  theDaughterName1,
112  theDaughterName2,
113  theDaughterName3),
114  parentmass(theParentMass),
115  theDaughterMasses(masses)
116 {
117  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
118 }
119 
120 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay(const G4String& theParentName,
121  G4double theParentMass,
122  G4double theBR,
123  G4int theNumberOfDaughters,
124  const G4String& theDaughterName1,
125  const G4String& theDaughterName2,
126  const G4String& theDaughterName3,
127  const G4String& theDaughterName4,
128  const G4double *masses) :
129  G4VDecayChannel("Phase Space",
130  theParentName,theBR,
131  theNumberOfDaughters,
132  theDaughterName1,
133  theDaughterName2,
134  theDaughterName3,
135  theDaughterName4),
136  parentmass(theParentMass),
137  theDaughterMasses(masses)
138 {
139  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
140 }
141 
142 G4GeneralPhaseSpaceDecay::~G4GeneralPhaseSpaceDecay()
143 {
144 }
145 
146 G4DecayProducts *G4GeneralPhaseSpaceDecay::DecayIt(G4double)
147 {
148  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
149  G4DecayProducts * products = NULL;
150 
151  CheckAndFillParent();
152  CheckAndFillDaughters();
153 
154  switch (numberOfDaughters){
155  case 0:
156  if (GetVerboseLevel()>0) {
157  G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
158  G4cout << " daughters not defined " <<G4endl;
159  }
160  break;
161  case 1:
162  products = OneBodyDecayIt();
163  break;
164  case 2:
165  products = TwoBodyDecayIt();
166  break;
167  case 3:
168  products = ThreeBodyDecayIt();
169  break;
170  default:
171  products = ManyBodyDecayIt();
172  break;
173  }
174  if ((products == NULL) && (GetVerboseLevel()>0)) {
175  G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
176  G4cout << *parent_name << " can not decay " << G4endl;
177  DumpInfo();
178  }
179  return products;
180 }
181 
182 G4DecayProducts *G4GeneralPhaseSpaceDecay::OneBodyDecayIt()
183 {
184  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::OneBodyDecayIt()"<<G4endl;
185 
186 // G4double daughtermass = daughters[0]->GetPDGMass();
187 
188  //create parent G4DynamicParticle at rest
189  G4ParticleMomentum dummy;
190  G4DynamicParticle * parentparticle = new G4DynamicParticle(G4MT_parent, dummy, 0.0);
191 
192  //create G4Decayproducts
193  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
194  delete parentparticle;
195 
196  //create daughter G4DynamicParticle at rest
197  G4DynamicParticle * daughterparticle = new G4DynamicParticle(G4MT_daughters[0], dummy, 0.0);
198  products->PushProducts(daughterparticle);
199 
200  if (GetVerboseLevel()>1)
201  {
202  G4cout << "G4GeneralPhaseSpaceDecay::OneBodyDecayIt ";
203  G4cout << " create decay products in rest frame " <<G4endl;
204  products->DumpInfo();
205  }
206  return products;
207 }
208 
209 G4DecayProducts *G4GeneralPhaseSpaceDecay::TwoBodyDecayIt()
210 {
211  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::TwoBodyDecayIt()"<<G4endl;
212 
213  //daughters'mass
214  G4double daughtermass[2];
215  G4double daughtermomentum;
216  if ( theDaughterMasses )
217  {
218  daughtermass[0]= *(theDaughterMasses);
219  daughtermass[1] = *(theDaughterMasses+1);
220  } else {
221  daughtermass[0] = G4MT_daughters[0]->GetPDGMass();
222  daughtermass[1] = G4MT_daughters[1]->GetPDGMass();
223  }
224 
225 // G4double sumofdaughtermass = daughtermass[0] + daughtermass[1];
226 
227  //create parent G4DynamicParticle at rest
228  G4ParticleMomentum dummy;
229  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
230 
231  //create G4Decayproducts @@GF why dummy parentparticle?
232  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
233  delete parentparticle;
234 
235  //calculate daughter momentum
236  daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
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  G4ParticleMomentum direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
241 
242  //create daughter G4DynamicParticle
243  G4double Etotal= std::sqrt(daughtermass[0]*daughtermass[0] + daughtermomentum*daughtermomentum);
244  G4DynamicParticle * daughterparticle = new G4DynamicParticle( G4MT_daughters[0],Etotal, direction*daughtermomentum);
245  products->PushProducts(daughterparticle);
246  Etotal= std::sqrt(daughtermass[1]*daughtermass[1] + daughtermomentum*daughtermomentum);
247  daughterparticle = new G4DynamicParticle( G4MT_daughters[1],Etotal, direction*(-1.0*daughtermomentum));
248  products->PushProducts(daughterparticle);
249 
250  if (GetVerboseLevel()>1)
251  {
252  G4cout << "G4GeneralPhaseSpaceDecay::TwoBodyDecayIt ";
253  G4cout << " create decay products in rest frame " <<G4endl;
254  products->DumpInfo();
255  }
256  return products;
257 }
258 
259 G4DecayProducts *G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt()
260 // algorism of this code is originally written in GDECA3 of GEANT3
261 {
262  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt()"<<G4endl;
263 
264  //daughters'mass
265  G4double daughtermass[3];
266  G4double sumofdaughtermass = 0.0;
267  for (G4int index=0; index<3; index++)
268  {
269  if ( theDaughterMasses )
270  {
271  daughtermass[index]= *(theDaughterMasses+index);
272  } else {
273  daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
274  }
275  sumofdaughtermass += daughtermass[index];
276  }
277 
278  //create parent G4DynamicParticle at rest
279  G4ParticleMomentum dummy;
280  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
281 
282  //create G4Decayproducts
283  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
284  delete parentparticle;
285 
286  //calculate daughter momentum
287  // Generate two
288  G4double rd1, rd2, rd;
289  G4double daughtermomentum[3];
290  G4double momentummax=0.0, momentumsum = 0.0;
291  G4double energy;
292  const G4int maxNumberOfLoops = 10000;
293  G4int loopCounter = 0;
294 
295  do
296  {
297  rd1 = G4UniformRand();
298  rd2 = G4UniformRand();
299  if (rd2 > rd1)
300  {
301  rd = rd1;
302  rd1 = rd2;
303  rd2 = rd;
304  }
305  momentummax = 0.0;
306  momentumsum = 0.0;
307  // daughter 0
308 
309  energy = rd2*(parentmass - sumofdaughtermass);
310  daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
311  if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0];
312  momentumsum += daughtermomentum[0];
313 
314  // daughter 1
315  energy = (1.-rd1)*(parentmass - sumofdaughtermass);
316  daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
317  if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1];
318  momentumsum += daughtermomentum[1];
319 
320  // daughter 2
321  energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
322  daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
323  if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2];
324  momentumsum += daughtermomentum[2];
325  } while ( ( momentummax > momentumsum - momentummax ) && /* Loop checking, 02.11.2015, A.Ribon */
326  ++loopCounter < maxNumberOfLoops );
327  if ( loopCounter >= maxNumberOfLoops ) {
328  G4ExceptionDescription ed;
329  ed << " Failed sampling after maxNumberOfLoops attempts : forced exit" << G4endl;
330  G4Exception( " G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt ", "HAD_PHASESPACE_001", FatalException, ed );
331  }
332 
333  // output message
334  if (GetVerboseLevel()>1) {
335  G4cout << " daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl;
336  G4cout << " daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl;
337  G4cout << " daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl;
338  G4cout << " momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl;
339  }
340 
341  //create daughter G4DynamicParticle
342  G4double costheta, sintheta, phi, sinphi, cosphi;
343  G4double costhetan, sinthetan, phin, sinphin, cosphin;
344  costheta = 2.*G4UniformRand()-1.0;
345  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
346  phi = twopi*G4UniformRand()*rad;
347  sinphi = std::sin(phi);
348  cosphi = std::cos(phi);
349  G4ParticleMomentum direction0(sintheta*cosphi,sintheta*sinphi,costheta);
350  G4double Etotal=std::sqrt( daughtermass[0]*daughtermass[0] + daughtermomentum[0]*daughtermomentum[0]);
351  G4DynamicParticle * daughterparticle
352  = new G4DynamicParticle( G4MT_daughters[0], Etotal, direction0*daughtermomentum[0]);
353  products->PushProducts(daughterparticle);
354 
355  costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
356  sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
357  phin = twopi*G4UniformRand()*rad;
358  sinphin = std::sin(phin);
359  cosphin = std::cos(phin);
360  G4ParticleMomentum direction2;
361  direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
362  direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
363  direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
364  Etotal=std::sqrt( daughtermass[2]*daughtermass[2] + daughtermomentum[2]*daughtermomentum[2]/direction2.mag2());
365  daughterparticle = new G4DynamicParticle( G4MT_daughters[2],Etotal, direction2*(daughtermomentum[2]/direction2.mag()));
366  products->PushProducts(daughterparticle);
367  G4ThreeVector mom=(direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0);
368  Etotal= std::sqrt( daughtermass[1]*daughtermass[1] + mom.mag2() );
369  daughterparticle =
370  new G4DynamicParticle(G4MT_daughters[1], Etotal, mom);
371  products->PushProducts(daughterparticle);
372 
373  if (GetVerboseLevel()>1) {
374  G4cout << "G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt ";
375  G4cout << " create decay products in rest frame " <<G4endl;
376  products->DumpInfo();
377  }
378  return products;
379 }
380 
381 G4DecayProducts *G4GeneralPhaseSpaceDecay::ManyBodyDecayIt()
382 // algorism of this code is originally written in FORTRAN by M.Asai
383 //*****************************************************************
384 // NBODY
385 // N-body phase space Monte-Carlo generator
386 // Makoto Asai
387 // Hiroshima Institute of Technology
388 // (asai@kekvax.kek.jp)
389 // Revised release : 19/Apr/1995
390 //
391 {
392  //return value
393  G4DecayProducts *products;
394 
395  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt()"<<G4endl;
396 
397  //daughters'mass
398  G4double *daughtermass = new G4double[numberOfDaughters];
399  G4double sumofdaughtermass = 0.0;
400  for (G4int index=0; index<numberOfDaughters; index++){
401  daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
402  sumofdaughtermass += daughtermass[index];
403  }
404 
405  //Calculate daughter momentum
406  G4double *daughtermomentum = new G4double[numberOfDaughters];
407  G4ParticleMomentum direction;
408  G4DynamicParticle **daughterparticle;
409  G4double *sm = new G4double[numberOfDaughters];
410  G4double tmas;
411  G4double weight = 1.0;
412  G4int numberOfTry = 0;
413  G4int index1;
414 
415  do {
416  //Generate rundom number in descending order
417  G4double temp;
418  G4double *rd = new G4double[numberOfDaughters];
419  rd[0] = 1.0;
420  for(index1 =1; index1 < numberOfDaughters -1; index1++)
421  rd[index1] = G4UniformRand();
422  rd[ numberOfDaughters -1] = 0.0;
423  for(index1 =1; index1 < numberOfDaughters -1; index1++) {
424  for(G4int index2 = index1+1; index2 < numberOfDaughters; index2++) {
425  if (rd[index1] < rd[index2]){
426  temp = rd[index1];
427  rd[index1] = rd[index2];
428  rd[index2] = temp;
429  }
430  }
431  }
432 
433  //calcurate virtual mass
434  tmas = parentmass - sumofdaughtermass;
435  temp = sumofdaughtermass;
436  for(index1 =0; index1 < numberOfDaughters; index1++) {
437  sm[index1] = rd[index1]*tmas + temp;
438  temp -= daughtermass[index1];
439  if (GetVerboseLevel()>1) {
440  G4cout << index1 << " rundom number:" << rd[index1];
441  G4cout << " virtual mass:" << sm[index1]/GeV << "[GeV/c/c]" <<G4endl;
442  }
443  }
444  delete [] rd;
445 
446  //Calculate daughter momentum
447  weight = 1.0;
448  index1 =numberOfDaughters-1;
449  daughtermomentum[index1]= Pmx( sm[index1-1],daughtermass[index1-1],sm[index1]);
450  if (GetVerboseLevel()>1) {
451  G4cout << " daughter " << index1 << ":" << *daughters_name[index1];
452  G4cout << " momentum:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
453  }
454  for(index1 =numberOfDaughters-2; index1>=0; index1--) {
455  // calculate
456  daughtermomentum[index1]= Pmx( sm[index1],daughtermass[index1], sm[index1 +1]);
457  if(daughtermomentum[index1] < 0.0) {
458  // !!! illegal momentum !!!
459  if (GetVerboseLevel()>0) {
460  G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt ";
461  G4cout << " can not calculate daughter momentum " <<G4endl;
462  G4cout << " parent:" << *parent_name;
463  G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
464  G4cout << " daughter " << index1 << ":" << *daughters_name[index1];
465  G4cout << " mass:" << daughtermass[index1]/GeV << "[GeV/c/c]" ;
466  G4cout << " mass:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
467  }
468  delete [] sm;
469  delete [] daughtermass;
470  delete [] daughtermomentum;
471  return NULL; // Error detection
472 
473  } else {
474  // calculate weight of this events
475  weight *= daughtermomentum[index1]/sm[index1];
476  if (GetVerboseLevel()>1) {
477  G4cout << " daughter " << index1 << ":" << *daughters_name[index1];
478  G4cout << " momentum:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
479  }
480  }
481  }
482  if (GetVerboseLevel()>1) {
483  G4cout << " weight: " << weight <<G4endl;
484  }
485 
486  // exit if number of Try exceeds 100
487  if (numberOfTry++ >100) {
488  if (GetVerboseLevel()>0) {
489  G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt: ";
490  G4cout << " can not determine Decay Kinematics " << G4endl;
491  }
492  delete [] sm;
493  delete [] daughtermass;
494  delete [] daughtermomentum;
495  return NULL; // Error detection
496  }
497  } while ( weight > G4UniformRand()); /* Loop checking, 02.11.2015, A.Ribon */
498  if (GetVerboseLevel()>1) {
499  G4cout << "Start calulation of daughters momentum vector "<<G4endl;
500  }
501 
502  G4double costheta, sintheta, phi;
503  G4double beta;
504  daughterparticle = new G4DynamicParticle*[numberOfDaughters];
505 
506  index1 = numberOfDaughters -2;
507  costheta = 2.*G4UniformRand()-1.0;
508  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
509  phi = twopi*G4UniformRand()*rad;
510  direction.setZ(costheta);
511  direction.setY(sintheta*std::sin(phi));
512  direction.setX(sintheta*std::cos(phi));
513  daughterparticle[index1] = new G4DynamicParticle( G4MT_daughters[index1], direction*daughtermomentum[index1] );
514  daughterparticle[index1+1] = new G4DynamicParticle( G4MT_daughters[index1+1], direction*(-1.0*daughtermomentum[index1]) );
515 
516  for (index1 = numberOfDaughters -3; index1 >= 0; index1--) {
517  //calculate momentum direction
518  costheta = 2.*G4UniformRand()-1.0;
519  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
520  phi = twopi*G4UniformRand()*rad;
521  direction.setZ(costheta);
522  direction.setY(sintheta*std::sin(phi));
523  direction.setX(sintheta*std::cos(phi));
524 
525  // boost already created particles
526  beta = daughtermomentum[index1];
527  beta /= std::sqrt( daughtermomentum[index1]*daughtermomentum[index1] + sm[index1+1]*sm[index1+1] );
528  for (G4int index2 = index1+1; index2<numberOfDaughters; index2++) {
529  G4LorentzVector p4;
530  // make G4LorentzVector for secondaries
531  p4 = daughterparticle[index2]->Get4Momentum();
532 
533  // boost secondaries to new frame
534  p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
535 
536  // change energy/momentum
537  daughterparticle[index2]->Set4Momentum(p4);
538  }
539  //create daughter G4DynamicParticle
540  daughterparticle[index1]= new G4DynamicParticle( G4MT_daughters[index1], direction*(-1.0*daughtermomentum[index1]));
541  }
542 
543  //create G4Decayproducts
544  G4DynamicParticle *parentparticle;
545  direction.setX(1.0); direction.setY(0.0); direction.setZ(0.0);
546  parentparticle = new G4DynamicParticle( G4MT_parent, direction, 0.0);
547  products = new G4DecayProducts(*parentparticle);
548  delete parentparticle;
549  for (index1 = 0; index1<numberOfDaughters; index1++) {
550  products->PushProducts(daughterparticle[index1]);
551  }
552  if (GetVerboseLevel()>1) {
553  G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt ";
554  G4cout << " create decay products in rest frame " << G4endl;
555  products->DumpInfo();
556  }
557 
558  delete [] daughterparticle;
559  delete [] daughtermomentum;
560  delete [] daughtermass;
561  delete [] sm;
562 
563  return products;
564 }
565 
566 
567 
568 
569