source: trunk/source/particles/management/src/G4PhaseSpaceDecayChannel.cc@ 1152

Last change on this file since 1152 was 992, checked in by garnier, 17 years ago

fichiers oublies
File size: 19.7 KB
Line 
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,v 1.13 2007/10/06 06:49:29 kurasige Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-02 $
29//
30//
31// ------------------------------------------------------------
32// GEANT 4 class header file
33//
34// History: first implementation, based on object model of
35// 27 July 1996 H.Kurashige
36// 20 Oct 1996 H.Kurashige
37// 30 May 1997 H.Kurashige
38// ------------------------------------------------------------
39
40#include "G4ParticleDefinition.hh"
41#include "G4DecayProducts.hh"
42#include "G4VDecayChannel.hh"
43#include "G4PhaseSpaceDecayChannel.hh"
44#include "Randomize.hh"
45#include "G4LorentzVector.hh"
46#include "G4LorentzRotation.hh"
47
48
49G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(G4int Verbose)
50 :G4VDecayChannel("Phase Space", Verbose)
51{
52
53}
54
55G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(
56 const G4String& theParentName,
57 G4double theBR,
58 G4int theNumberOfDaughters,
59 const G4String& theDaughterName1,
60 const G4String& theDaughterName2,
61 const G4String& theDaughterName3,
62 const G4String& theDaughterName4 )
63 :G4VDecayChannel("Phase Space",
64 theParentName,theBR,
65 theNumberOfDaughters,
66 theDaughterName1,
67 theDaughterName2,
68 theDaughterName3,
69 theDaughterName4)
70{
71
72}
73
74G4PhaseSpaceDecayChannel::~G4PhaseSpaceDecayChannel()
75{
76}
77
78G4DecayProducts *G4PhaseSpaceDecayChannel::DecayIt(G4double parentMass)
79{
80#ifdef G4VERBOSE
81 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
82#endif
83
84 G4DecayProducts * products = 0;
85
86 if (parent == 0) FillParent();
87 if (daughters == 0) FillDaughters();
88
89 if (parentMass >0.0) parent_mass = parentMass;
90 switch (numberOfDaughters){
91 case 0:
92#ifdef G4VERBOSE
93 if (GetVerboseLevel()>0) {
94 G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
95 G4cout << " daughters not defined " <<G4endl;
96 }
97#endif
98 break;
99 case 1:
100 products = OneBodyDecayIt();
101 break;
102 case 2:
103 products = TwoBodyDecayIt();
104 break;
105 case 3:
106 products = ThreeBodyDecayIt();
107 break;
108 default:
109 products = ManyBodyDecayIt();
110 break;
111 }
112#ifdef G4VERBOSE
113 if ((products == 0) && (GetVerboseLevel()>0)) {
114 G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
115 G4cout << *parent_name << " can not decay " << G4endl;
116 DumpInfo();
117 }
118#endif
119 return products;
120}
121
122G4DecayProducts *G4PhaseSpaceDecayChannel::OneBodyDecayIt()
123{
124#ifdef G4VERBOSE
125 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt()"<<G4endl;
126#endif
127
128 //create parent G4DynamicParticle at rest
129 G4ThreeVector dummy;
130 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
131 //create G4Decayproducts
132 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
133 delete parentparticle;
134
135 //create daughter G4DynamicParticle at rest
136 G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0], dummy, 0.0);
137 products->PushProducts(daughterparticle);
138
139#ifdef G4VERBOSE
140 if (GetVerboseLevel()>1) {
141 G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt ";
142 G4cout << " create decay products in rest frame " <<G4endl;
143 products->DumpInfo();
144 }
145#endif
146 return products;
147}
148
149G4DecayProducts *G4PhaseSpaceDecayChannel::TwoBodyDecayIt()
150{
151#ifdef G4VERBOSE
152 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt()"<<G4endl;
153#endif
154 // parent mass
155 G4double parentmass = parent_mass;
156
157 //daughters'mass
158 G4double daughtermass[2];
159 G4double daughtermomentum;
160 daughtermass[0] = daughters_mass[0];
161 daughtermass[1] = daughters_mass[1];
162
163 //create parent G4DynamicParticle at rest
164 G4ThreeVector dummy;
165 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
166 //create G4Decayproducts
167 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
168 delete parentparticle;
169
170 //calculate daughter momentum
171 daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
172 if (daughtermomentum <0.0) {
173#ifdef G4VERBOSE
174 if (GetVerboseLevel()>0) {
175 G4cerr << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt "
176 << "sum of daughter mass is larger than parent mass" << G4endl;
177 G4cerr << "parent :" << parent->GetParticleName() << " " << parent_mass/GeV << G4endl;
178 G4cerr << "daughter 1 :" << daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
179 G4cerr << "daughter 2:" << daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
180 }
181#endif
182 G4Exception("G4PhaseSpaceDecayChannel::TwoBodyDecayIt",
183 "can not create decay products", JustWarning,
184 "sum of daughter mass is larger than parent mass");
185 return products;
186 }
187
188 G4double costheta = 2.*G4UniformRand()-1.0;
189 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
190 G4double phi = twopi*G4UniformRand()*rad;
191 G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
192
193 //create daughter G4DynamicParticle
194 G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0], direction*daughtermomentum);
195 products->PushProducts(daughterparticle);
196 daughterparticle = new G4DynamicParticle( daughters[1], direction*(-1.0*daughtermomentum));
197 products->PushProducts(daughterparticle);
198
199#ifdef G4VERBOSE
200 if (GetVerboseLevel()>1) {
201 G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt ";
202 G4cout << " create decay products in rest frame " <<G4endl;
203 products->DumpInfo();
204 }
205#endif
206 return products;
207}
208
209G4DecayProducts *G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()
210// algorism of this code is originally written in GDECA3 of GEANT3
211{
212#ifdef G4VERBOSE
213 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()"<<G4endl;
214#endif
215 // parent mass
216 G4double parentmass = parent_mass;
217 //daughters'mass
218 G4double daughtermass[3];
219 G4double sumofdaughtermass = 0.0;
220 for (G4int index=0; index<3; index++){
221 daughtermass[index] = daughters_mass[index];
222 sumofdaughtermass += daughtermass[index];
223 }
224
225 //create parent G4DynamicParticle at rest
226 G4ThreeVector dummy;
227 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
228
229
230 //create G4Decayproducts
231 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
232 delete parentparticle;
233
234 if (sumofdaughtermass >parentmass) {
235#ifdef G4VERBOSE
236 if (GetVerboseLevel()>0) {
237 G4cerr << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt "
238 << "sum of daughter mass is larger than parent mass" << G4endl;
239 G4cerr << "parent :" << parent->GetParticleName() << " " << parent_mass/GeV << G4endl;
240 G4cerr << "daughter 1 :" << daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
241 G4cerr << "daughter 2:" << daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
242 G4cerr << "daughter 3:" << daughters[2]->GetParticleName() << " " << daughtermass[2]/GeV << G4endl;
243 }
244#endif
245 G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
246 "can not create decay products", JustWarning,
247 "sum of daughter mass is larger than parent mass");
248 return products;
249 }
250
251
252
253 //calculate daughter momentum
254 // Generate two
255 G4double rd1, rd2, rd;
256 G4double daughtermomentum[3];
257 G4double momentummax=0.0, momentumsum = 0.0;
258 G4double energy;
259
260 do {
261 rd1 = G4UniformRand();
262 rd2 = G4UniformRand();
263 if (rd2 > rd1) {
264 rd = rd1;
265 rd1 = rd2;
266 rd2 = rd;
267 }
268 momentummax = 0.0;
269 momentumsum = 0.0;
270 // daughter 0
271 energy = rd2*(parentmass - sumofdaughtermass);
272 daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
273 if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0];
274 momentumsum += daughtermomentum[0];
275 // daughter 1
276 energy = (1.-rd1)*(parentmass - sumofdaughtermass);
277 daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
278 if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1];
279 momentumsum += daughtermomentum[1];
280 // daughter 2
281 energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
282 daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
283 if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2];
284 momentumsum += daughtermomentum[2];
285 } while (momentummax > momentumsum - momentummax );
286 // output message
287#ifdef G4VERBOSE
288 if (GetVerboseLevel()>1) {
289 G4cout << " daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl;
290 G4cout << " daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl;
291 G4cout << " daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl;
292 G4cout << " momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl;
293 }
294#endif
295
296 //create daughter G4DynamicParticle
297 G4double costheta, sintheta, phi, sinphi, cosphi;
298 G4double costhetan, sinthetan, phin, sinphin, cosphin;
299 costheta = 2.*G4UniformRand()-1.0;
300 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
301 phi = twopi*G4UniformRand()*rad;
302 sinphi = std::sin(phi);
303 cosphi = std::cos(phi);
304 G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
305 G4DynamicParticle * daughterparticle
306 = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
307 products->PushProducts(daughterparticle);
308
309 costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
310 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
311 phin = twopi*G4UniformRand()*rad;
312 sinphin = std::sin(phin);
313 cosphin = std::cos(phin);
314 G4ThreeVector direction2;
315 direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
316 direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
317 direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
318 daughterparticle = new G4DynamicParticle( daughters[2], direction2*(daughtermomentum[2]/direction2.mag()));
319 products->PushProducts(daughterparticle);
320
321 daughterparticle =
322 new G4DynamicParticle(
323 daughters[1],
324 (direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0)
325 );
326 products->PushProducts(daughterparticle);
327
328#ifdef G4VERBOSE
329 if (GetVerboseLevel()>1) {
330 G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt ";
331 G4cout << " create decay products in rest frame " <<G4endl;
332 products->DumpInfo();
333 }
334#endif
335 return products;
336}
337
338G4DecayProducts *G4PhaseSpaceDecayChannel::ManyBodyDecayIt()
339// algorithm of this code is originally written in FORTRAN by M.Asai
340//******************************************************************
341// NBODY
342// N-body phase space Monte-Carlo generator
343// Makoto Asai
344// Hiroshima Institute of Technology
345// (asai@kekvax.kek.jp)
346// Revised release : 19/Apr/1995
347//
348{
349 G4int index, index2;
350
351 //return value
352 G4DecayProducts *products;
353
354#ifdef G4VERBOSE
355 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()"<<G4endl;
356#endif
357 // parent mass
358 G4double parentmass = parent_mass;
359 //daughters'mass
360 G4double *daughtermass = new G4double[numberOfDaughters];
361 G4double sumofdaughtermass = 0.0;
362 for (index=0; index<numberOfDaughters; index++){
363 daughtermass[index] = daughters_mass[index];
364 sumofdaughtermass += daughtermass[index];
365 }
366
367 //Calculate daughter momentum
368 G4double *daughtermomentum = new G4double[numberOfDaughters];
369 G4ThreeVector direction;
370 G4DynamicParticle **daughterparticle;
371 G4double *sm = new G4double[numberOfDaughters];
372 G4double tmas;
373 G4double weight = 1.0;
374 G4int numberOfTry = 0;
375
376 do {
377 //Generate rundom number in descending order
378 G4double temp;
379 G4double *rd = new G4double[numberOfDaughters];
380 rd[0] = 1.0;
381 for(index =1; index < numberOfDaughters -1; index++) rd[index] = G4UniformRand();
382 rd[ numberOfDaughters -1] = 0.0;
383 for(index =1; index < numberOfDaughters -1; index++) {
384 for(index2 = index+1; index2 < numberOfDaughters; index2++) {
385 if (rd[index] < rd[index2]){
386 temp = rd[index];
387 rd[index] = rd[index2];
388 rd[index2] = temp;
389 }
390 }
391 }
392
393 //calcurate virtual mass
394 tmas = parentmass - sumofdaughtermass;
395 temp = sumofdaughtermass;
396 for(index =0; index < numberOfDaughters; index++) {
397 sm[index] = rd[index]*tmas + temp;
398 temp -= daughtermass[index];
399 if (GetVerboseLevel()>1) {
400 G4cout << index << " rundom number:" << rd[index];
401 G4cout << " virtual mass:" << sm[index]/GeV << "[GeV/c/c]" <<G4endl;
402 }
403 }
404 delete [] rd;
405
406 //Calculate daughter momentum
407 weight = 1.0;
408 index =numberOfDaughters-1;
409 daughtermomentum[index]= Pmx( sm[index-1],daughtermass[index-1], sm[index]);
410#ifdef G4VERBOSE
411 if (GetVerboseLevel()>1) {
412 G4cout << " daughter " << index << ":" << *daughters_name[index];
413 G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
414 }
415#endif
416 for(index =numberOfDaughters-2; index>=0; index--) {
417 // calculate
418 daughtermomentum[index]= Pmx( sm[index],daughtermass[index], sm[index +1]);
419 if(daughtermomentum[index] < 0.0) {
420 // !!! illegal momentum !!!
421#ifdef G4VERBOSE
422 if (GetVerboseLevel()>0) {
423 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
424 G4cout << " can not calculate daughter momentum " <<G4endl;
425 G4cout << " parent:" << *parent_name;
426 G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
427 G4cout << " daughter " << index << ":" << *daughters_name[index];
428 G4cout << " mass:" << daughtermass[index]/GeV << "[GeV/c/c]" ;
429 G4cout << " mass:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
430 }
431#endif
432 delete [] sm;
433 delete [] daughtermass;
434 delete [] daughtermomentum;
435
436 G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt",
437 "can not create decay products", JustWarning,
438 "sum of daughter mass is larger than parent mass");
439
440 return 0; // Error detection
441
442 } else {
443 // calculate weight of this events
444 weight *= daughtermomentum[index]/sm[index];
445#ifdef G4VERBOSE
446 if (GetVerboseLevel()>1) {
447 G4cout << " daughter " << index << ":" << *daughters_name[index];
448 G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
449 }
450#endif
451 }
452 }
453
454#ifdef G4VERBOSE
455 if (GetVerboseLevel()>1) {
456 G4cout << " weight: " << weight <<G4endl;
457 }
458#endif
459
460 // exit if number of Try exceeds 100
461 if (numberOfTry++ >100) {
462#ifdef G4VERBOSE
463 if (GetVerboseLevel()>0) {
464 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ";
465 G4cout << " can not determine Decay Kinematics " << G4endl;
466 G4cout << "parent : " << *parent_name << G4endl;
467 G4cout << "daughters : ";
468 for(index=0; index<numberOfDaughters; index++) {
469 G4cout << *daughters_name[index] << " , ";
470 }
471 G4cout << G4endl;
472 }
473#endif
474 G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ",
475 " Cannot decay ", JustWarning,
476 " Decay Kinematics cannot be calculated ");
477
478 delete [] sm;
479 delete [] daughtermass;
480 delete [] daughtermomentum;
481 return 0; // Error detection
482 }
483 } while ( weight > G4UniformRand());
484
485#ifdef G4VERBOSE
486 if (GetVerboseLevel()>1) {
487 G4cout << "Start calulation of daughters momentum vector "<<G4endl;
488 }
489#endif
490
491 G4double costheta, sintheta, phi;
492 G4double beta;
493 daughterparticle = new G4DynamicParticle*[numberOfDaughters];
494
495 index = numberOfDaughters -2;
496 costheta = 2.*G4UniformRand()-1.0;
497 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
498 phi = twopi*G4UniformRand()*rad;
499 direction.setZ(costheta);
500 direction.setY(sintheta*std::sin(phi));
501 direction.setX(sintheta*std::cos(phi));
502 daughterparticle[index] = new G4DynamicParticle( daughters[index], direction*daughtermomentum[index] );
503 daughterparticle[index+1] = new G4DynamicParticle( daughters[index+1], direction*(-1.0*daughtermomentum[index]) );
504
505 for (index = numberOfDaughters -3; index >= 0; index--) {
506 //calculate momentum direction
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
514 // boost already created particles
515 beta = daughtermomentum[index];
516 beta /= std::sqrt( daughtermomentum[index]*daughtermomentum[index] + sm[index+1]*sm[index+1] );
517 for (G4int index2 = index+1; index2<numberOfDaughters; index2++) {
518 G4LorentzVector p4;
519 // make G4LorentzVector for secondaries
520 p4 = daughterparticle[index2]->Get4Momentum();
521
522 // boost secondaries to new frame
523 p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
524
525 // change energy/momentum
526 daughterparticle[index2]->Set4Momentum(p4);
527 }
528 //create daughter G4DynamicParticle
529 daughterparticle[index]= new G4DynamicParticle( daughters[index], direction*(-1.0*daughtermomentum[index]));
530 }
531
532 //create G4Decayproducts
533 G4DynamicParticle *parentparticle;
534 direction.setX(1.0); direction.setY(0.0); direction.setZ(0.0);
535 parentparticle = new G4DynamicParticle( parent, direction, 0.0);
536 products = new G4DecayProducts(*parentparticle);
537 delete parentparticle;
538 for (index = 0; index<numberOfDaughters; index++) {
539 products->PushProducts(daughterparticle[index]);
540 }
541
542#ifdef G4VERBOSE
543 if (GetVerboseLevel()>1) {
544 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
545 G4cout << " create decay products in rest frame " << G4endl;
546 products->DumpInfo();
547 }
548#endif
549
550 delete [] daughterparticle;
551 delete [] daughtermomentum;
552 delete [] daughtermass;
553 delete [] sm;
554
555 return products;
556}
557
558
559G4double G4PhaseSpaceDecayChannel::Pmx(G4double e, G4double p1, G4double p2)
560{
561 // calcurate momentum of daughter particles in two-body decay
562 G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e);
563 if (ppp>0) return std::sqrt(ppp);
564 else return -1.;
565}
566
567
568
Note: See TracBrowser for help on using the repository browser.