| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 // GEANT 4 class header file
29 //
30 // History: first implementation, A. Feli
31 //
32 // Note: this class is a generalization o
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::G4GeneralPhaseSpaceD
49 G4VDecayChannel("Pha
50 parentmass(0.), theD
51 {
52 if (GetVerboseLevel()>1) G4cout << "G4Genera
53 }
54
55 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceD
56 G4double
57 G4int
58 const G4Strin
59 const G4Strin
60 const G4Strin
61 G4VDecayCha
62 theParentName,theBR,
63 theNumberOfDaughters,
64 theDaughterName1,
65 theDaughterName2,
66 theDaughterName3),
67 theDaughterMasses(0)
68 {
69 if (GetVerboseLevel()>1) G4cout << "G4Genera
70
71 // Set the parent particle (resonance) mas
72 if (G4MT_parent != NULL)
73 {
74 parentmass = G4MT_parent->GetPDGMass();
75 } else {
76 parentmass=0.;
77 }
78
79 }
80
81 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceD
82
83 G4double
84 G4int
85 const G4Strin
86 const G4Strin
87 const G4Strin
88 G4VDecayCha
89 theParentName,theBR,
90 theNumberOfDaughters,
91 theDaughterName1,
92 theDaughterName2,
93 theDaughterName3),
94 parentmass(theParentMass),
95 theDaughterMasses(0)
96 {
97 if (GetVerboseLevel()>1) G4cout << "G4Genera
98 }
99
100 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceD
101
102 G4double
103 G4int
104 const G4Strin
105 const G4Strin
106 const G4Strin
107 const G4double *masses) :
108 G4VDecayCha
109 theParentName,theBR,
110 theNumberOfDaughters,
111 theDaughterName1,
112 theDaughterName2,
113 theDaughterName3),
114 parentmass(theParentMass),
115 theDaughterMasses(masses)
116 {
117 if (GetVerboseLevel()>1) G4cout << "G4Genera
118 }
119
120 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceD
121
122 G4double
123 G4int
124 const G4Strin
125 const G4Strin
126 const G4Strin
127 const G4Strin
128 const G4double *masses) :
129 G4VDecayCha
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 << "G4Genera
140 }
141
142 G4GeneralPhaseSpaceDecay::~G4GeneralPhaseSpace
143 {
144 }
145
146 G4DecayProducts *G4GeneralPhaseSpaceDecay::Dec
147 {
148 if (GetVerboseLevel()>1) G4cout << "G4Genera
149 G4DecayProducts * products = NULL;
150
151 CheckAndFillParent();
152 CheckAndFillDaughters();
153
154 switch (numberOfDaughters){
155 case 0:
156 if (GetVerboseLevel()>0) {
157 G4cout << "G4GeneralPhaseSpaceDecay::Dec
158 G4cout << " daughters not defined " <<G4
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()
175 G4cout << "G4GeneralPhaseSpaceDecay::Decay
176 G4cout << *parent_name << " can not decay
177 DumpInfo();
178 }
179 return products;
180 }
181
182 G4DecayProducts *G4GeneralPhaseSpaceDecay::One
183 {
184 if (GetVerboseLevel()>1) G4cout << "G4Genera
185
186 // G4double daughtermass = daughters[0]->GetP
187
188 //create parent G4DynamicParticle at rest
189 G4ParticleMomentum dummy;
190 G4DynamicParticle * parentparticle = new G4D
191
192 //create G4Decayproducts
193 G4DecayProducts *products = new G4DecayProdu
194 delete parentparticle;
195
196 //create daughter G4DynamicParticle at rest
197 G4DynamicParticle * daughterparticle = new G
198 products->PushProducts(daughterparticle);
199
200 if (GetVerboseLevel()>1)
201 {
202 G4cout << "G4GeneralPhaseSpaceDecay::OneB
203 G4cout << " create decay products in res
204 products->DumpInfo();
205 }
206 return products;
207 }
208
209 G4DecayProducts *G4GeneralPhaseSpaceDecay::Two
210 {
211 if (GetVerboseLevel()>1) G4cout << "G4Genera
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]->GetP
222 daughtermass[1] = G4MT_daughters[1]->GetP
223 }
224
225 // G4double sumofdaughtermass = daughtermass
226
227 //create parent G4DynamicParticle at rest
228 G4ParticleMomentum dummy;
229 G4DynamicParticle * parentparticle = new G4D
230
231 //create G4Decayproducts @@GF why dummy par
232 G4DecayProducts *products = new G4DecayProdu
233 delete parentparticle;
234
235 //calculate daughter momentum
236 daughtermomentum = Pmx(parentmass,daughterma
237 G4double costheta = 2.*G4UniformRand()-1.0;
238 G4double sintheta = std::sqrt((1.0 - costhet
239 G4double phi = twopi*G4UniformRand()*rad;
240 G4ParticleMomentum direction(sintheta*std::c
241
242 //create daughter G4DynamicParticle
243 G4double Etotal= std::sqrt(daughtermass[0]*d
244 G4DynamicParticle * daughterparticle = new G
245 products->PushProducts(daughterparticle);
246 Etotal= std::sqrt(daughtermass[1]*daughterma
247 daughterparticle = new G4DynamicParticle( G4
248 products->PushProducts(daughterparticle);
249
250 if (GetVerboseLevel()>1)
251 {
252 G4cout << "G4GeneralPhaseSpaceDecay::TwoB
253 G4cout << " create decay products in res
254 products->DumpInfo();
255 }
256 return products;
257 }
258
259 G4DecayProducts *G4GeneralPhaseSpaceDecay::Thr
260 // algorism of this code is originally written
261 {
262 if (GetVerboseLevel()>1) G4cout << "G4Genera
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]= *(theDaughterMas
272 } else {
273 daughtermass[index] = G4MT_daughters[
274 }
275 sumofdaughtermass += daughtermass[index];
276 }
277
278 //create parent G4DynamicParticle at rest
279 G4ParticleMomentum dummy;
280 G4DynamicParticle * parentparticle = new G4D
281
282 //create G4Decayproducts
283 G4DecayProducts *products = new G4DecayProdu
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 - sumofdaughterm
310 daughtermomentum[0] = std::sqrt(energy*en
311 if ( daughtermomentum[0] >momentummax )mo
312 momentumsum += daughtermomentum[0];
313
314 // daughter 1
315 energy = (1.-rd1)*(parentmass - sumofdaug
316 daughtermomentum[1] = std::sqrt(energy*en
317 if ( daughtermomentum[1] >momentummax )mo
318 momentumsum += daughtermomentum[1];
319
320 // daughter 2
321 energy = (rd1-rd2)*(parentmass - sumofdau
322 daughtermomentum[2] = std::sqrt(energy*en
323 if ( daughtermomentum[2] >momentummax )mo
324 momentumsum += daughtermomentum[2];
325 } while ( ( momentummax > momentumsum - m
326 ++loopCounter < maxNumberOfLoops
327 if ( loopCounter >= maxNumberOfLoops ) {
328 G4ExceptionDescription ed;
329 ed << " Failed sampling after maxNumberO
330 G4Exception( " G4GeneralPhaseSpaceDecay:
331 }
332
333 // output message
334 if (GetVerboseLevel()>1) {
335 G4cout << " daughter 0:" << daughtermo
336 G4cout << " daughter 1:" << daughtermo
337 G4cout << " daughter 2:" << daughtermo
338 G4cout << " momentum sum:" << momentumsu
339 }
340
341 //create daughter G4DynamicParticle
342 G4double costheta, sintheta, phi, sinphi, co
343 G4double costhetan, sinthetan, phin, sinphin
344 costheta = 2.*G4UniformRand()-1.0;
345 sintheta = std::sqrt((1.0-costheta)*(1.0+cos
346 phi = twopi*G4UniformRand()*rad;
347 sinphi = std::sin(phi);
348 cosphi = std::cos(phi);
349 G4ParticleMomentum direction0(sintheta*cosph
350 G4double Etotal=std::sqrt( daughtermass[0]*d
351 G4DynamicParticle * daughterparticle
352 = new G4DynamicParticle( G4MT_daughte
353 products->PushProducts(daughterparticle);
354
355 costhetan = (daughtermomentum[1]*daughtermom
356 sinthetan = std::sqrt((1.0-costhetan)*(1.0+c
357 phin = twopi*G4UniformRand()*rad;
358 sinphin = std::sin(phin);
359 cosphin = std::cos(phin);
360 G4ParticleMomentum direction2;
361 direction2.setX( sinthetan*cosphin*costheta*
362 direction2.setY( sinthetan*cosphin*costheta*
363 direction2.setZ( -sinthetan*cosphin*sintheta
364 Etotal=std::sqrt( daughtermass[2]*daughterma
365 daughterparticle = new G4DynamicParticle( G4
366 products->PushProducts(daughterparticle);
367 G4ThreeVector mom=(direction0*daughtermoment
368 Etotal= std::sqrt( daughtermass[1]*daughterm
369 daughterparticle =
370 new G4DynamicParticle(G4MT_daughters[1]
371 products->PushProducts(daughterparticle);
372
373 if (GetVerboseLevel()>1) {
374 G4cout << "G4GeneralPhaseSpaceDecay::Thre
375 G4cout << " create decay products in res
376 products->DumpInfo();
377 }
378 return products;
379 }
380
381 G4DecayProducts *G4GeneralPhaseSpaceDecay::Man
382 // algorism of this code is originally written
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 << "G4Genera
396
397 //daughters'mass
398 G4double *daughtermass = new G4double[number
399 G4double sumofdaughtermass = 0.0;
400 for (G4int index=0; index<numberOfDaughters;
401 daughtermass[index] = G4MT_daughters[index
402 sumofdaughtermass += daughtermass[index];
403 }
404
405 //Calculate daughter momentum
406 G4double *daughtermomentum = new G4double[nu
407 G4ParticleMomentum direction;
408 G4DynamicParticle **daughterparticle;
409 G4double *sm = new G4double[numberOfDaughter
410 G4double tmas;
411 G4double weight = 1.0;
412 G4int numberOfTry = 0;
413 G4int index1;
414
415 do {
416 //Generate rundom number in descending ord
417 G4double temp;
418 G4double *rd = new G4double[numberOfDaught
419 rd[0] = 1.0;
420 for(index1 =1; index1 < numberOfDaughters
421 rd[index1] = G4UniformRand();
422 rd[ numberOfDaughters -1] = 0.0;
423 for(index1 =1; index1 < numberOfDaughters
424 for(G4int index2 = index1+1; index2 < nu
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;
437 sm[index1] = rd[index1]*tmas + temp;
438 temp -= daughtermass[index1];
439 if (GetVerboseLevel()>1) {
440 G4cout << index1 << " rundom number:"
441 G4cout << " virtual mass:" << sm[ind
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
450 if (GetVerboseLevel()>1) {
451 G4cout << " daughter " << index1 <<
452 G4cout << " momentum:" << daughtermoment
453 }
454 for(index1 =numberOfDaughters-2; index1>=0
455 // calculate
456 daughtermomentum[index1]= Pmx( sm[index1
457 if(daughtermomentum[index1] < 0.0) {
458 // !!! illegal momentum !!!
459 if (GetVerboseLevel()>0) {
460 G4cout << "G4GeneralPhaseSpaceDecay:
461 G4cout << " can not calculate da
462 G4cout << " parent:" << *parent_
463 G4cout << " mass:" << parentmass/GeV
464 G4cout << " daughter " << index1
465 G4cout << " mass:" << daughtermass[i
466 G4cout << " mass:" << daughtermoment
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
476 if (GetVerboseLevel()>1) {
477 G4cout << " daughter " << index1
478 G4cout << " momentum:" << daughtermo
479 }
480 }
481 }
482 if (GetVerboseLevel()>1) {
483 G4cout << " weight: " << weight <<G4e
484 }
485
486 // exit if number of Try exceeds 100
487 if (numberOfTry++ >100) {
488 if (GetVerboseLevel()>0) {
489 G4cout << "G4GeneralPhaseSpaceDecay::M
490 G4cout << " can not determine Decay Kinemati
491 }
492 delete [] sm;
493 delete [] daughtermass;
494 delete [] daughtermomentum;
495 return NULL; // Error detection
496 }
497 } while ( weight > G4UniformRand()); /* Loo
498 if (GetVerboseLevel()>1) {
499 G4cout << "Start calculation of daughter
500 }
501
502 G4double costheta, sintheta, phi;
503 G4double beta;
504 daughterparticle = new G4DynamicParticle*[nu
505
506 index1 = numberOfDaughters -2;
507 costheta = 2.*G4UniformRand()-1.0;
508 sintheta = std::sqrt((1.0-costheta)*(1.0+cos
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 G4DynamicPart
514 daughterparticle[index1+1] = new G4DynamicPa
515
516 for (index1 = numberOfDaughters -3; index1
517 //calculate momentum direction
518 costheta = 2.*G4UniformRand()-1.0;
519 sintheta = std::sqrt((1.0-costheta)*(1.0+c
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
528 for (G4int index2 = index1+1; index2<numbe
529 G4LorentzVector p4;
530 // make G4LorentzVector for secondaries
531 p4 = daughterparticle[index2]->Get4Momen
532
533 // boost secondaries to new frame
534 p4.boost( direction.x()*beta, direction.
535
536 // change energy/momentum
537 daughterparticle[index2]->Set4Momentum(p
538 }
539 //create daughter G4DynamicParticle
540 daughterparticle[index1]= new G4DynamicPar
541 }
542
543 //create G4Decayproducts
544 G4DynamicParticle *parentparticle;
545 direction.setX(1.0); direction.setY(0.0); d
546 parentparticle = new G4DynamicParticle( G4MT
547 products = new G4DecayProducts(*parentpartic
548 delete parentparticle;
549 for (index1 = 0; index1<numberOfDaughters; i
550 products->PushProducts(daughterparticle[in
551 }
552 if (GetVerboseLevel()>1) {
553 G4cout << "G4GeneralPhaseSpaceDecay::ManyB
554 G4cout << " create decay products in rest
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