source: trunk/source/processes/hadronic/models/util/src/G4GeneralPhaseSpaceDecay.cc@ 1201

Last change on this file since 1201 was 819, checked in by garnier, 17 years ago

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