| Geant4 Cross Reference |
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 calculation 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
570