Main interface to the de-excitation code for hyper-nuclei.
92{
93
96
98
99mult10:
101
102 varntp->clear();
103
104 if (nucleusS > 0)
105 nucleusS = 0;
106
107 G4int NbLam0 = std::abs(nucleusS);
108
109 Ainit = -1 * nucleusA;
110 Zinit = -1 * nucleusZ;
111 Sinit = -1 * nucleusS;
112
115 G4int ZFP1 = 0, AFP1 = 0, AFPIMF = 0, ZFPIMF = 0, ZFP2 = 0, AFP2 = 0, SFP1 = 0, SFP2 = 0, SFPIMF = 0;
116 G4double vx_eva = 0.0, vy_eva = 0.0, vz_eva = 0.0;
117 G4double VX_PREF = 0., VY_PREF = 0., VZ_PREF = 00, VP1X, VP1Y, VP1Z, VXOUT, VYOUT, VZOUT, V_CM[3], VFP1_CM[3],
118 VFP2_CM[3], VIMF_CM[3], VX2OUT, VY2OUT, VZ2OUT;
119 G4double zf = 0.0, af = 0.0, mtota = 0.0, tkeimf = 0.0, jprf0 = 0.;
120 G4int ff = 0, afpnew = 0, zfpnew = 0, aprfp = 0, zprfp = 0, IOUNSTABLE = 0, ILOOP = 0, IEV_TAB = 0,
121 IEV_TAB_TEMP = 0;
122 G4int fimf = 0, INMIN = 0, INMAX = 0;
124 G4int inum = eventnumber;
126 opt->optimfallowed = 1;
127
128 if (fiss->zt > 56)
129 {
130 fiss->ifis = 1;
131 }
132 else
133 {
134 fiss->ifis = 0;
135 }
136
137 if (NbLam0 > 0)
138 {
139 opt->nblan0 = NbLam0;
140 }
141
146
151
152 gammaemission = 0;
153 G4double T_init = 0., T_diff = 0., a_tilda = 0., a_tilda_BU = 0., EE_diff = 0., EINCL = 0., A_FINAL = 0.,
154 Z_FINAL = 0., E_FINAL = 0.;
155
156 G4double A_diff = 0., ASLOPE1, ASLOPE2, A_ACC, ABU_SLOPE, ABU_SUM = 0., AMEM = 0., ZMEM = 0., EMEM = 0., JMEM = 0.,
157 PX_BU_SUM = 0.0, PY_BU_SUM = 0.0, PZ_BU_SUM = 0.0, ETOT_SUM = 0., P_BU_SUM = 0., ZBU_SUM = 0.,
158 Z_Breakup_sum = 0., A_Breakup, Z_Breakup, N_Breakup, G_SYMM, CZ, Sigma_Z, Z_Breakup_Mean, ZTEMP = 0.,
159 ATEMP = 0.;
160
161 G4double ETOT_PRF = 0.0, PXPRFP = 0., PYPRFP = 0., PZPRFP = 0., PPRFP = 0., VX1_BU = 0., VY1_BU = 0., VZ1_BU = 0.,
162 VBU2 = 0., GAMMA_REL = 1.0, Eexc_BU_SUM = 0., VX_BU_SUM = 0., VY_BU_SUM = 0., VZ_BU_SUM = 0.,
163 E_tot_BU = 0., EKIN_BU = 0., ZIMFBU = 0., AIMFBU = 0., ZFFBU = 0., AFFBU = 0., AFBU = 0., ZFBU = 0.,
164 EEBU = 0., TKEIMFBU = 0., vx_evabu = 0., vy_evabu = 0., vz_evabu = 0., Bvalue_BU = 0., P_BU = 0.,
165 ETOT_BU = 1., PX_BU = 0., PY_BU = 0., PZ_BU = 0., VX2_BU = 0., VY2_BU = 0., VZ2_BU = 0.;
166
167 G4int ABU_DIFF, ZBU_DIFF, NBU_DIFF;
168 G4int INEWLOOP = 0, ILOOPBU = 0;
169
173
174 if (nucleusA < 1)
175 {
176 std::cout << "Error - Remnant with a mass number A below 1." << std::endl;
177
178 return;
179 }
180
181 for (
G4int j = 0; j < 3; j++)
182 {
183 V_CM[j] = 0.;
184 VFP1_CM[j] = 0.;
185 VFP2_CM[j] = 0.;
186 VIMF_CM[j] = 0.;
187 }
188
190 {
191 for (
G4int I2 = 0; I2 < 12; I2++)
192 BU_TAB[I1][I2] = 0.0;
193 for (
G4int I2 = 0; I2 < 6; I2++)
194 {
195 BU_TAB_TEMP[I1][I2] = 0.0;
196 BU_TAB_TEMP1[I1][I2] = 0.0;
197 EV_TAB_TEMP[I1][I2] = 0.0;
198 EV_TAB[I1][I2] = 0.0;
199 EV_TAB_SSC[I1][I2] = 0.0;
200 EV_TEMP[I1][I2] = 0.0;
201 }
202 }
203
205 if (idebug == 1)
206 {
207 zprf = 81.;
208 aprf = 201.;
209
210 ee = 100.0;
211 jprf = 10.;
212 zf = 0.;
213 af = 0.;
214 mtota = 0.;
215 ff = 1;
216 inttype = 0;
217
218 }
219
222 EINCL = ee;
223
224
225
226
227
228 G4double pincl = std::sqrt(pxrem * pxrem + pyrem * pyrem + pzrem * pzrem);
229
230 G4double ETOT_incl = std::sqrt(pincl * pincl + (AAINCL * amu) * (AAINCL * amu));
231 G4double VX_incl =
C * pxrem / ETOT_incl;
232 G4double VY_incl =
C * pyrem / ETOT_incl;
233 G4double VZ_incl =
C * pzrem / ETOT_incl;
234
235
240 IEV_TAB = 0;
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263 if (T_freeze_out_in >= 0.0)
264 {
265 T_freeze_out = T_freeze_out_in;
266 }
267 else
268 {
269 T_freeze_out =
max(9.33 * std::exp(-0.00282 * AAINCL), 5.5);
270
271
272
273 }
274
275 a_tilda = ald->av * aprf + ald->as * std::pow(aprf, 2.0 / 3.0) + ald->ak * std::pow(aprf, 1.0 / 3.0);
276
277 T_init = std::sqrt(EINCL / a_tilda);
278
279 T_diff = T_init - T_freeze_out;
280
281 if (T_diff > 0.1 && zprf > 2. && (aprf - zprf) > 0.)
282 {
283
284
285 varntp->kfis = 10;
286
287 for (
G4int i = 0; i < 5; i++)
288 {
289 EE_diff = EINCL - a_tilda * T_freeze_out * T_freeze_out;
290
291
292
293
294
295
296 A_diff =
dint(EE_diff / (8.0 * 5.0 / T_freeze_out));
297
298 if (A_diff > AAINCL)
299 A_diff = AAINCL;
300
301 A_FINAL = AAINCL - A_diff;
302
303 a_tilda =
304 ald->av * A_FINAL + ald->as * std::pow(A_FINAL, 2.0 / 3.0) + ald->ak * std::pow(A_FINAL, 1.0 / 3.0);
305 E_FINAL = a_tilda * T_freeze_out * T_freeze_out;
306
307 if (A_FINAL < 4.0)
308 {
309 EE_diff = EINCL - E_FINAL;
310 A_FINAL = 1.0;
311 Z_FINAL = 1.0;
312 E_FINAL = 0.0;
313 goto mul4325;
314 }
315 }
316 mul4325:
317
318
319
320
321 Z_FINAL =
dint(zprf * A_FINAL / (aprf));
322
323 if (E_FINAL < 0.0)
324 E_FINAL = 0.0;
325
326 aprf = A_FINAL;
327 zprf = Z_FINAL;
328 ee = E_FINAL;
329
330 A_diff = AAINCL - aprf;
331
332
333 if (A_diff <= 1.0)
334 {
335 aprf = AAINCL;
336 zprf = ZAINCL;
337 ee = EINCL;
338 IMULTIFR = 0;
339 goto mult7777;
340 }
341 else if (A_diff > 1.0)
342 {
343
344 A_ACC = 0.0;
345
346
347 ASLOPE1 = -2.400;
348 ASLOPE2 = -1.200;
349
350 a_tilda = ald->av * AAINCL + ald->as * std::pow(AAINCL, 2.0 / 3.0) + ald->ak * std::pow(AAINCL, 1.0 / 3.0);
351
352 E_FINAL = a_tilda * T_freeze_out * T_freeze_out;
353
354 ABU_SLOPE = (ASLOPE1 - ASLOPE2) / 4.0 * (E_FINAL / AAINCL) + ASLOPE1 - (ASLOPE1 - ASLOPE2) * 7.0 / 4.0;
355
356
357
358
359
360
361
362
363
364
365
366 if (ABU_SLOPE > -1.01)
367 ABU_SLOPE = -1.01;
368
369 I_Breakup = 0;
370 Z_Breakup_sum = Z_FINAL;
371 ABU_SUM = 0.0;
372 ZBU_SUM = 0.0;
373
374 for (
G4int i = 0; i < 100; i++)
375 {
376 IS = 0;
377 mult4326:
379
380 IS = IS + 1;
381 if (IS > 100)
382 {
383 std::cout << "WARNING: IPOWERLIMHAZ CALLED MORE THAN 100 TIMES WHEN "
384 "CALCULATING A_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED: "
385 << A_Breakup << std::endl;
386 goto mult10;
387 }
388
389 if (A_Breakup > AAINCL)
390 goto mult4326;
391
392 if (A_Breakup <= 0.0)
393 {
394 std::cout << "A_BREAKUP <= 0 " << std::endl;
395 goto mult10;
396 }
397
398 A_ACC = A_ACC + A_Breakup;
399
400 if (A_ACC <= A_diff)
401 {
402
403 Z_Breakup_Mean =
dint(A_Breakup * ZAINCL / AAINCL);
404
405 Z_Breakup_sum = Z_Breakup_sum + Z_Breakup_Mean;
406
407
408 G_SYMM = 34.2281 - 5.14037 * E_FINAL / AAINCL;
409 if (E_FINAL / AAINCL < 2.0)
410 G_SYMM = 25.0;
411 if (E_FINAL / AAINCL > 4.0)
412 G_SYMM = 15.0;
413
414
415
416 G_SYMM = 25.0;
417 CZ = 2.0 * G_SYMM * 4.0 / A_Breakup;
418
419
421 Sigma_Z = std::sqrt(T_freeze_out / CZ);
422
423 IS = 0;
424 mult4333:
426 IS = IS + 1;
427
428 if (IS > 100)
429 {
430 std::cout << "WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN "
431 "CALCULATING Z_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE "
432 "DICED: "
433 << A_Breakup << " " << Z_Breakup << std::endl;
434 goto mult10;
435 }
436
437 if (Z_Breakup < 0.0)
438 goto mult4333;
439 if ((A_Breakup - Z_Breakup) < 0.0)
440 goto mult4333;
441 if ((A_Breakup - Z_Breakup) == 0.0 && Z_Breakup != 1.0)
442 goto mult4333;
443
444 if (Z_Breakup >= ZAINCL)
445 {
446 IIS = IIS + 1;
447 if (IIS > 10)
448 {
449 std::cout << "Z_BREAKUP RESAMPLED MORE THAN 10 TIMES; EVENT WILL "
450 "BE RESAMPLED AGAIN "
451 << std::endl;
452 goto mult10;
453 }
454 goto mult4333;
455 }
456
457
459
460 if (Z_Breakup > 2.0)
461 {
462 if (
idnint(A_Breakup - Z_Breakup) < INMIN ||
idnint(A_Breakup - Z_Breakup) > (INMAX + 5))
463 {
464
465
466 goto mult4343;
467 }
468 }
469
470 mult4343:
471
472
473
474
475 N_Breakup = A_Breakup - Z_Breakup;
476 BU_TAB[I_Breakup][0] =
dint(Z_Breakup);
477 BU_TAB[I_Breakup][1] =
dint(A_Breakup);
478 ABU_SUM = ABU_SUM + BU_TAB[i][1];
479 ZBU_SUM = ZBU_SUM + BU_TAB[i][0];
480
481
482 BU_TAB[I_Breakup][3] = 0.0;
483 I_Breakup = I_Breakup + 1;
484 IMULTBU = IMULTBU + 1;
485 }
486 else
487 {
488
489
490
491
492
493
494 goto mult4327;
495 }
496 }
497
498
499 }
500 mult4327:
501 IMULTIFR = 1;
502
503
504 ABU_DIFF =
idnint(ABU_SUM + aprf - AAINCL);
505 ZBU_DIFF =
idnint(ZBU_SUM + zprf - ZAINCL);
506 NBU_DIFF =
idnint((ABU_SUM - ZBU_SUM) + (aprf - zprf) - (AAINCL - ZAINCL));
507
508 if (IMULTBU > 200)
509 std::cout << "WARNING - MORE THAN 200 BU " << IMULTBU << std::endl;
510
511 if (IMULTBU < 1)
512 std::cout << "WARNING - LESS THAN 1 BU " << IMULTBU << std::endl;
513
514
516 for (
G4int i = 0; i < IMULTBU; i++)
517 IMEM_BU[i] = 0;
518
519 while (NBU_DIFF != 0 && ZBU_DIFF != 0)
520 {
521
522
523
524 IS = 0;
525 mult5555:
527 IPROBA = IPROBA + 1;
528 IS = IS + 1;
529 if (IS > 100)
530 {
531 std::cout << "WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING "
532 "N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED."
533 << std::endl;
534 goto mult10;
535 }
537 if (IMEM_BU[IEL] == 1)
538 goto mult5555;
539 if (!(IEL < 200))
540 std::cout << "5555:" << IEL << RHAZ << IMULTBU << std::endl;
541 if (IEL < 0)
542 std::cout << "5555:" << IEL << RHAZ << IMULTBU << std::endl;
543 if (IEL <= IMULTBU)
544 {
545 N_Breakup =
dint(BU_TAB[IEL][1] - BU_TAB[IEL][0] -
DSIGN(1.0,
G4double(NBU_DIFF)));
546 }
547 else if (IEL > IMULTBU)
548 {
550 }
551 if (N_Breakup < 0.0)
552 {
553 IMEM_BU[IEL] = 1;
554 goto mult5555;
555 }
556 if (IEL <= IMULTBU)
557 {
559 }
560 else if (IEL > IMULTBU)
561 {
563 }
564 if (ZTEMP < 0.0)
565 {
566 IMEM_BU[IEL] = 1;
567 goto mult5555;
568 }
569 if (ZTEMP < 1.0 && N_Breakup < 1.0)
570 {
571 IMEM_BU[IEL] = 1;
572 goto mult5555;
573 }
574
575
576
577
578
579
580
581 if (IEL <= IMULTBU)
582 {
583 BU_TAB[IEL][0] =
dint(ZTEMP);
584 BU_TAB[IEL][1] =
dint(ZTEMP + N_Breakup);
585 }
586 else if (IEL > IMULTBU)
587 {
589 aprf =
dint(ZTEMP + N_Breakup);
590 }
591 NBU_DIFF = NBU_DIFF -
ISIGN(1, NBU_DIFF);
592 ZBU_DIFF = ZBU_DIFF -
ISIGN(1, ZBU_DIFF);
593 }
594
595 IPROBA = 0;
596 for (
G4int i = 0; i < IMULTBU; i++)
597 IMEM_BU[i] = 0;
598
599 if (NBU_DIFF != 0 && ZBU_DIFF == 0)
600 {
601 while (NBU_DIFF > 0 || NBU_DIFF < 0)
602 {
603 IS = 0;
604 mult5556:
606 IS = IS + 1;
607 if (IS > 100)
608 {
609 std::cout << "WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING "
610 "N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED."
611 << std::endl;
612 goto mult10;
613 }
615 if (IMEM_BU[IEL] == 1)
616 goto mult5556;
617
618 if (IPROBA > IMULTBU + 1 && NBU_DIFF > 0)
619 {
620 std::cout << "###',IPROBA,IMULTBU,NBU_DIFF,ZBU_DIFF,T_freeze_out" << std::endl;
621 IPROBA = IPROBA + 1;
622 if (IEL <= IMULTBU)
623 {
624 BU_TAB[IEL][1] =
dint(BU_TAB[IEL][1] -
G4double(NBU_DIFF));
625 }
626 else
627 {
628 if (IEL > IMULTBU)
630 }
631 goto mult5432;
632 }
633 if (!(IEL < 200))
634 std::cout << "5556:" << IEL << RHAZ << IMULTBU << std::endl;
635 if (IEL < 0)
636 std::cout << "5556:" << IEL << RHAZ << IMULTBU << std::endl;
637 if (IEL <= IMULTBU)
638 {
639 N_Breakup =
dint(BU_TAB[IEL][1] - BU_TAB[IEL][0] -
DSIGN(1.0,
G4double(NBU_DIFF)));
640 }
641 else if (IEL > IMULTBU)
642 {
644 }
645 if (N_Breakup < 0.0)
646 {
647 IMEM_BU[IEL] = 1;
648 goto mult5556;
649 }
650 if (IEL <= IMULTBU)
651 {
652 ATEMP =
dint(BU_TAB[IEL][0] + N_Breakup);
653 }
654 else if (IEL > IMULTBU)
655 {
656 ATEMP =
dint(zprf + N_Breakup);
657 }
658 if ((ATEMP - N_Breakup) < 1.0 && N_Breakup < 1.0)
659 {
660 IMEM_BU[IEL] = 1;
661 goto mult5556;
662 }
663
664
665
666
667
668 if (IEL <= IMULTBU)
669 BU_TAB[IEL][1] =
dint(BU_TAB[IEL][0] + N_Breakup);
670 else if (IEL > IMULTBU)
671 aprf =
dint(zprf + N_Breakup);
672
673 NBU_DIFF = NBU_DIFF -
ISIGN(1, NBU_DIFF);
674 }
675
676 IPROBA = 0;
677 for (
G4int i = 0; i < IMULTBU; i++)
678 IMEM_BU[i] = 0;
679 }
680 else
681 {
682 if (ZBU_DIFF != 0 && NBU_DIFF == 0)
683 {
684 while (ZBU_DIFF > 0 || ZBU_DIFF < 0)
685 {
686 IS = 0;
687 mult5557:
689 IS = IS + 1;
690 if (IS > 100)
691 {
692 std::cout << "WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING "
693 "N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED."
694 << std::endl;
695 goto mult10;
696 }
698 if (IMEM_BU[IEL] == 1)
699 goto mult5557;
700
701 if (IPROBA > IMULTBU + 1 && ZBU_DIFF > 0)
702 {
703 std::cout << "###',IPROBA,IMULTBU,NBU_DIFF,ZBU_DIFF,T_freeze_out" << std::endl;
704 IPROBA = IPROBA + 1;
705 if (IEL <= IMULTBU)
706 {
707 N_Breakup =
dint(BU_TAB[IEL][1] - BU_TAB[IEL][0]);
708 BU_TAB[IEL][0] =
dint(BU_TAB[IEL][0] -
G4double(ZBU_DIFF));
709 BU_TAB[IEL][1] =
dint(BU_TAB[IEL][0] + N_Breakup);
710 }
711 else
712 {
713 if (IEL > IMULTBU)
714 {
715 N_Breakup = aprf - zprf;
717 aprf =
dint(zprf + N_Breakup);
718 }
719 }
720 goto mult5432;
721 }
722 if (!(IEL < 200))
723 std::cout << "5557:" << IEL << RHAZ << IMULTBU << std::endl;
724 if (IEL < 0)
725 std::cout << "5557:" << IEL << RHAZ << IMULTBU << std::endl;
726 if (IEL <= IMULTBU)
727 {
728 N_Breakup =
dint(BU_TAB[IEL][1] - BU_TAB[IEL][0]);
730 }
731 else if (IEL > IMULTBU)
732 {
733 N_Breakup =
dint(aprf - zprf);
735 }
736 ATEMP =
dint(ZTEMP + N_Breakup);
737 if (ZTEMP < 0.0)
738 {
739 IMEM_BU[IEL] = 1;
740 goto mult5557;
741 }
742 if ((ATEMP - ZTEMP) < 0.0)
743 {
744 IMEM_BU[IEL] = 1;
745 goto mult5557;
746 }
747 if ((ATEMP - ZTEMP) < 1.0 && ZTEMP < 1.0)
748 {
749 IMEM_BU[IEL] = 1;
750 goto mult5557;
751 }
752 if (IEL <= IMULTBU)
753 {
754 BU_TAB[IEL][0] =
dint(ZTEMP);
755 BU_TAB[IEL][1] =
dint(ZTEMP + N_Breakup);
756 }
757 else
758 {
759 if (IEL > IMULTBU)
760 {
762 aprf =
dint(ZTEMP + N_Breakup);
763 }
764 }
765 ZBU_DIFF = ZBU_DIFF -
ISIGN(1, ZBU_DIFF);
766 }
767 }
768 }
769
770 mult5432:
771
772
773 ZMEM = 0.0;
774
775 for (
G4int i = 0; i < IMULTBU; i++)
776 {
777
778
779
780
781
782 if (BU_TAB[i][0] > 2.0)
783 {
784 a_tilda_BU = ald->av * BU_TAB[i][1] + ald->as * std::pow(BU_TAB[i][1], 2.0 / 3.0) +
785 ald->ak * std::pow(BU_TAB[i][1], 1.0 / 3.0);
786 BU_TAB[i][2] = a_tilda_BU * T_freeze_out * T_freeze_out;
787 }
788 else
789 {
790 BU_TAB[i][2] = 0.0;
791 }
792
793 if (BU_TAB[i][0] > ZMEM)
794 {
795 IMEM = i;
796 ZMEM = BU_TAB[i][0];
797 AMEM = BU_TAB[i][1];
798 EMEM = BU_TAB[i][2];
799 JMEM = BU_TAB[i][3];
800 }
801 }
802
803 if (zprf < ZMEM)
804 {
805 BU_TAB[IMEM][0] = zprf;
806 BU_TAB[IMEM][1] = aprf;
807 BU_TAB[IMEM][2] = ee;
808 BU_TAB[IMEM][3] = jprf;
809 zprf = ZMEM;
810 aprf = AMEM;
813 ee = EMEM;
814 jprf = JMEM;
815 }
816
817
818 ABU_SUM = aprf;
819 ZBU_SUM = zprf;
820 for (
G4int i = 0; i < IMULTBU; i++)
821 {
822 ABU_SUM = ABU_SUM + BU_TAB[i][1];
823 ZBU_SUM = ZBU_SUM + BU_TAB[i][0];
824 }
825 ABU_DIFF =
idnint(ABU_SUM - AAINCL);
826 ZBU_DIFF =
idnint(ZBU_SUM - ZAINCL);
827
828 if (ABU_DIFF != 0 || ZBU_DIFF != 0)
829 std::cout << "Problem of mass in BU " << ABU_DIFF << " " << ZBU_DIFF << std::endl;
830 PX_BU_SUM = 0.0;
831 PY_BU_SUM = 0.0;
832 PZ_BU_SUM = 0.0;
833
834
835
836 AMOMENT(AAINCL, aprf, 1, &PXPRFP, &PYPRFP, &PZPRFP);
837 PPRFP = std::sqrt(PXPRFP * PXPRFP + PYPRFP * PYPRFP + PZPRFP * PZPRFP);
838
839
840 ETOT_PRF = std::sqrt(PPRFP * PPRFP + (aprf * amu) * (aprf * amu));
841 VX_PREF =
C * PXPRFP / ETOT_PRF;
842 VY_PREF =
C * PYPRFP / ETOT_PRF;
843 VZ_PREF =
C * PZPRFP / ETOT_PRF;
844
845
846 tke_bu(zprf, aprf, ZAINCL, AAINCL, &VX1_BU, &VY1_BU, &VZ1_BU);
847
848
849
850
851
852
853 lorentz_boost(VX1_BU, VY1_BU, VZ1_BU, VX_PREF, VY_PREF, VZ_PREF, &VXOUT, &VYOUT, &VZOUT);
854
855 VX_PREF = VXOUT;
856 VY_PREF = VYOUT;
857 VZ_PREF = VZOUT;
858
859
860 VBU2 = VX_PREF * VX_PREF + VY_PREF * VY_PREF + VZ_PREF * VZ_PREF;
861 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C * C));
862 ETOT_PRF = aprf * amu / GAMMA_REL;
863 PXPRFP = ETOT_PRF * VX_PREF /
C;
864 PYPRFP = ETOT_PRF * VY_PREF /
C;
865 PZPRFP = ETOT_PRF * VZ_PREF /
C;
866
867
868
869
870
871
872 PX_BU_SUM = PXPRFP;
873 PY_BU_SUM = PYPRFP;
874 PZ_BU_SUM = PZPRFP;
875
876 Eexc_BU_SUM = ee;
878
879 for (I_Breakup = 0; I_Breakup < IMULTBU; I_Breakup++)
880 {
881
882 Bvalue_BU = Bvalue_BU +
eflmac(
idnint(BU_TAB[I_Breakup][1]),
idnint(BU_TAB[I_Breakup][0]), 1, 0);
883 Eexc_BU_SUM = Eexc_BU_SUM + BU_TAB[I_Breakup][2];
884
885 AMOMENT(AAINCL, BU_TAB[I_Breakup][1], 1, &PX_BU, &PY_BU, &PZ_BU);
886 P_BU = std::sqrt(PX_BU * PX_BU + PY_BU * PY_BU + PZ_BU * PZ_BU);
887
888
889 ETOT_BU = std::sqrt(P_BU * P_BU + (BU_TAB[I_Breakup][1] * amu) * (BU_TAB[I_Breakup][1] * amu));
890 BU_TAB[I_Breakup][4] =
C * PX_BU / ETOT_BU;
891 BU_TAB[I_Breakup][5] =
C * PY_BU / ETOT_BU;
892 BU_TAB[I_Breakup][6] =
C * PZ_BU / ETOT_BU;
893
894 tke_bu(BU_TAB[I_Breakup][0], BU_TAB[I_Breakup][1], ZAINCL, AAINCL, &VX2_BU, &VY2_BU, &VZ2_BU);
895
896
897
898
899
900
902 VY2_BU,
903 VZ2_BU,
904 BU_TAB[I_Breakup][4],
905 BU_TAB[I_Breakup][5],
906 BU_TAB[I_Breakup][6],
907 &VXOUT,
908 &VYOUT,
909 &VZOUT);
910
911 BU_TAB[I_Breakup][4] = VXOUT;
912 BU_TAB[I_Breakup][5] = VYOUT;
913 BU_TAB[I_Breakup][6] = VZOUT;
914
915
916 VBU2 = BU_TAB[I_Breakup][4] * BU_TAB[I_Breakup][4] + BU_TAB[I_Breakup][5] * BU_TAB[I_Breakup][5] +
917 BU_TAB[I_Breakup][6] * BU_TAB[I_Breakup][6];
918 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C * C));
919 ETOT_BU = BU_TAB[I_Breakup][1] * amu / GAMMA_REL;
920 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
921 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
922 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
923
924 PX_BU_SUM = PX_BU_SUM + PX_BU;
925 PY_BU_SUM = PY_BU_SUM + PY_BU;
926 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
927
928 }
929
930
931 P_BU_SUM = std::sqrt(PX_BU_SUM * PX_BU_SUM + PY_BU_SUM * PY_BU_SUM + PZ_BU_SUM * PZ_BU_SUM);
932
933
934 ETOT_SUM = std::sqrt(P_BU_SUM * P_BU_SUM + (AAINCL * amu) * (AAINCL * amu));
935
936 VX_BU_SUM =
C * PX_BU_SUM / ETOT_SUM;
937 VY_BU_SUM =
C * PY_BU_SUM / ETOT_SUM;
938 VZ_BU_SUM =
C * PZ_BU_SUM / ETOT_SUM;
939
940
941
942
943
944
945 lorentz_boost(-VX_BU_SUM, -VY_BU_SUM, -VZ_BU_SUM, VX_PREF, VY_PREF, VZ_PREF, &VXOUT, &VYOUT, &VZOUT);
946
947 VX_PREF = VXOUT;
948 VY_PREF = VYOUT;
949 VZ_PREF = VZOUT;
950
951 VBU2 = VX_PREF * VX_PREF + VY_PREF * VY_PREF + VZ_PREF * VZ_PREF;
952 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C * C));
953 ETOT_PRF = aprf * amu / GAMMA_REL;
954 PXPRFP = ETOT_PRF * VX_PREF /
C;
955 PYPRFP = ETOT_PRF * VY_PREF /
C;
956 PZPRFP = ETOT_PRF * VZ_PREF /
C;
957
958 PX_BU_SUM = 0.0;
959 PY_BU_SUM = 0.0;
960 PZ_BU_SUM = 0.0;
961
962 PX_BU_SUM = PXPRFP;
963 PY_BU_SUM = PYPRFP;
964 PZ_BU_SUM = PZPRFP;
965 E_tot_BU = ETOT_PRF;
966
967 EKIN_BU = aprf * amu / GAMMA_REL - aprf * amu;
968
969 for (I_Breakup = 0; I_Breakup < IMULTBU; I_Breakup++)
970 {
971
972
973
974
975
977 -VY_BU_SUM,
978 -VZ_BU_SUM,
979 BU_TAB[I_Breakup][4],
980 BU_TAB[I_Breakup][5],
981 BU_TAB[I_Breakup][6],
982 &VXOUT,
983 &VYOUT,
984 &VZOUT);
985
986 BU_TAB[I_Breakup][4] = VXOUT;
987 BU_TAB[I_Breakup][5] = VYOUT;
988 BU_TAB[I_Breakup][6] = VZOUT;
989
990 VBU2 = BU_TAB[I_Breakup][4] * BU_TAB[I_Breakup][4] + BU_TAB[I_Breakup][5] * BU_TAB[I_Breakup][5] +
991 BU_TAB[I_Breakup][6] * BU_TAB[I_Breakup][6];
992 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C * C));
993
994 ETOT_BU = BU_TAB[I_Breakup][1] * amu / GAMMA_REL;
995
996 EKIN_BU = EKIN_BU + BU_TAB[I_Breakup][1] * amu / GAMMA_REL - BU_TAB[I_Breakup][1] * amu;
997
998 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
999 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
1000 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
1001 E_tot_BU = E_tot_BU + ETOT_BU;
1002
1003 PX_BU_SUM = PX_BU_SUM + PX_BU;
1004 PY_BU_SUM = PY_BU_SUM + PY_BU;
1005 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
1006 }
1007
1008 if (std::abs(PX_BU_SUM) > 10. || std::abs(PY_BU_SUM) > 10. || std::abs(PZ_BU_SUM) > 10.)
1009 {
1010
1011
1012 P_BU_SUM = std::sqrt(PX_BU_SUM * PX_BU_SUM + PY_BU_SUM * PY_BU_SUM + PZ_BU_SUM * PZ_BU_SUM);
1013
1014
1015 ETOT_SUM = std::sqrt(P_BU_SUM * P_BU_SUM + (AAINCL * amu) * (AAINCL * amu));
1016
1017 VX_BU_SUM =
C * PX_BU_SUM / ETOT_SUM;
1018 VY_BU_SUM =
C * PY_BU_SUM / ETOT_SUM;
1019 VZ_BU_SUM =
C * PZ_BU_SUM / ETOT_SUM;
1020
1021
1022
1023
1024
1025
1026 lorentz_boost(-VX_BU_SUM, -VY_BU_SUM, -VZ_BU_SUM, VX_PREF, VY_PREF, VZ_PREF, &VXOUT, &VYOUT, &VZOUT);
1027
1028 VX_PREF = VXOUT;
1029 VY_PREF = VYOUT;
1030 VZ_PREF = VZOUT;
1031
1032 VBU2 = VX_PREF * VX_PREF + VY_PREF * VY_PREF + VZ_PREF * VZ_PREF;
1033 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C * C));
1034 ETOT_PRF = aprf * amu / GAMMA_REL;
1035 PXPRFP = ETOT_PRF * VX_PREF /
C;
1036 PYPRFP = ETOT_PRF * VY_PREF /
C;
1037 PZPRFP = ETOT_PRF * VZ_PREF /
C;
1038
1039 PX_BU_SUM = 0.0;
1040 PY_BU_SUM = 0.0;
1041 PZ_BU_SUM = 0.0;
1042
1043 PX_BU_SUM = PXPRFP;
1044 PY_BU_SUM = PYPRFP;
1045 PZ_BU_SUM = PZPRFP;
1046 E_tot_BU = ETOT_PRF;
1047
1048 EKIN_BU = aprf * amu / GAMMA_REL - aprf * amu;
1049
1050 for (I_Breakup = 0; I_Breakup < IMULTBU; I_Breakup++)
1051 {
1052
1053
1054
1055
1056
1058 -VY_BU_SUM,
1059 -VZ_BU_SUM,
1060 BU_TAB[I_Breakup][4],
1061 BU_TAB[I_Breakup][5],
1062 BU_TAB[I_Breakup][6],
1063 &VXOUT,
1064 &VYOUT,
1065 &VZOUT);
1066
1067 BU_TAB[I_Breakup][4] = VXOUT;
1068 BU_TAB[I_Breakup][5] = VYOUT;
1069 BU_TAB[I_Breakup][6] = VZOUT;
1070
1071 VBU2 = BU_TAB[I_Breakup][4] * BU_TAB[I_Breakup][4] + BU_TAB[I_Breakup][5] * BU_TAB[I_Breakup][5] +
1072 BU_TAB[I_Breakup][6] * BU_TAB[I_Breakup][6];
1073 GAMMA_REL = std::sqrt(1.0 - VBU2 / (C * C));
1074
1075 ETOT_BU = BU_TAB[I_Breakup][1] * amu / GAMMA_REL;
1076
1077 EKIN_BU = EKIN_BU + BU_TAB[I_Breakup][1] * amu / GAMMA_REL - BU_TAB[I_Breakup][1] * amu;
1078
1079 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
1080 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
1081 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
1082 E_tot_BU = E_tot_BU + ETOT_BU;
1083
1084 PX_BU_SUM = PX_BU_SUM + PX_BU;
1085 PY_BU_SUM = PY_BU_SUM + PY_BU;
1086 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
1087 }
1088 }
1089
1090
1091
1092
1093
1094 INEWLOOP = 0;
1095 for (
G4int i = 0; i < IMULTBU; i++)
1096 {
1097 if (BU_TAB[i][0] < 3.0 || BU_TAB[i][0] == BU_TAB[i][1])
1098 {
1101 &afpnew,
1102 &zfpnew,
1103 IOUNSTABLE,
1104 BU_TAB[i][4],
1105 BU_TAB[i][5],
1106 BU_TAB[i][6],
1107 &VP1X,
1108 &VP1Y,
1109 &VP1Z,
1110 BU_TAB_TEMP,
1111 &ILOOP);
1112
1113 if (IOUNSTABLE > 0)
1114 {
1115
1118 BU_TAB[i][4] = VP1X;
1119 BU_TAB[i][5] = VP1Y;
1120 BU_TAB[i][6] = VP1Z;
1121
1122
1123 for (int IJ = 0; IJ < ILOOP; IJ++)
1124 {
1125 BU_TAB[IMULTBU + INEWLOOP + IJ][0] = BU_TAB_TEMP[IJ][0];
1126 BU_TAB[IMULTBU + INEWLOOP + IJ][1] = BU_TAB_TEMP[IJ][1];
1127 BU_TAB[IMULTBU + INEWLOOP + IJ][4] = BU_TAB_TEMP[IJ][2];
1128 BU_TAB[IMULTBU + INEWLOOP + IJ][5] = BU_TAB_TEMP[IJ][3];
1129 BU_TAB[IMULTBU + INEWLOOP + IJ][6] = BU_TAB_TEMP[IJ][4];
1130 BU_TAB[IMULTBU + INEWLOOP + IJ][2] = 0.0;
1131 BU_TAB[IMULTBU + INEWLOOP + IJ][3] = 0.0;
1132 }
1133
1134 INEWLOOP = INEWLOOP + ILOOP;
1135
1136 }
1137 }
1138 }
1139
1140
1141 IMULTBU = IMULTBU + INEWLOOP;
1142
1143 opt->optimfallowed = 1;
1144 fiss->ifis = 0;
1145 gammaemission = 0;
1146 ILOOPBU = 0;
1147
1148
1152 for (
G4int i = 0; i < IMULTBU; i++)
1153 sumN = sumN + BU_TAB[i][1] - BU_TAB[i][0];
1154
1155 for (
G4int i = 0; i < IMULTBU; i++)
1156 {
1157 problamb[i] = (BU_TAB[i][1] - BU_TAB[i][0]) / sumN;
1158 }
1160 Nblamb =
new G4int[IMULTBU];
1161 for (
G4int i = 0; i < IMULTBU; i++)
1162 Nblamb[i] = 0;
1163 for (
G4int j = 0; j < NbLam0;)
1164 {
1165 G4double probtotal = (aprf - zprf) / sumN;
1167
1168 if (ran <= probtotal)
1169 {
1170 NbLamprf++;
1171 goto directlamb0;
1172 }
1173 for (
G4int i = 0; i < IMULTBU; i++)
1174 {
1175
1176 if (probtotal < ran && ran <= probtotal + problamb[i])
1177 {
1178 Nblamb[i] = Nblamb[i] + 1;
1179 goto directlamb0;
1180 }
1181 probtotal = probtotal + problamb[i];
1182 }
1183 directlamb0:
1184 j++;
1185 }
1186
1187 for (
G4int i = 0; i < IMULTBU; i++)
1188 {
1189 EEBU = BU_TAB[i][2];
1190 BU_TAB[i][10] = BU_TAB[i][6];
1192 if (BU_TAB[i][0] > 2.0)
1193 {
1194 G4int nbl = Nblamb[i];
1196 BU_TAB[i][1],
1197 &EEBU,
1198 0.0,
1199 &ZFBU,
1200 &AFBU,
1201 &mtota,
1202 &vz_evabu,
1203 &vx_evabu,
1204 &vy_evabu,
1205 &ff,
1206 &fimf,
1207 &ZIMFBU,
1208 &AIMFBU,
1209 &TKEIMFBU,
1210 &jprfbu,
1211 &inttype,
1212 &inum,
1213 EV_TEMP,
1214 &IEV_TAB_TEMP,
1215 &nbl);
1216
1217 Nblamb[i] = nbl;
1218 BU_TAB[i][9] = jprfbu;
1219
1220
1221
1222 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++)
1223 {
1224 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1225 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1226 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1227
1228
1229
1230
1231
1233 BU_TAB[i][5],
1234 BU_TAB[i][6],
1235 EV_TEMP[IJ][2],
1236 EV_TEMP[IJ][3],
1237 EV_TEMP[IJ][4],
1238 &VXOUT,
1239 &VYOUT,
1240 &VZOUT);
1241 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1242 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1243 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1244 }
1245 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1246
1247
1248
1249
1250
1251
1252
1254 vx_evabu, vy_evabu, vz_evabu, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VXOUT, &VYOUT, &VZOUT);
1255 BU_TAB[i][4] = VXOUT;
1256 BU_TAB[i][5] = VYOUT;
1257 BU_TAB[i][6] = VZOUT;
1258
1259 if (fimf == 0)
1260 {
1261 BU_TAB[i][7] =
dint(ZFBU);
1262 BU_TAB[i][8] =
dint(AFBU);
1263 BU_TAB[i][11] = nbl;
1264 }
1265
1266 if (fimf == 1)
1267 {
1268
1269
1270
1271
1274 opt->optimfallowed = 0;
1275 fiss->ifis = 0;
1276
1277 G4double EkinR1 = TKEIMFBU * AIMFBU / (AFBU + AIMFBU);
1278 G4double EkinR2 = TKEIMFBU * AFBU / (AFBU + AIMFBU);
1279 G4double V1 = std::sqrt(EkinR1 / AFBU) * 1.3887;
1280 G4double V2 = std::sqrt(EkinR2 / AIMFBU) * 1.3887;
1282 G4double VPERP1 = std::sqrt(
V1 *
V1 - VZ1_IMF * VZ1_IMF);
1284 G4double VX1_IMF = VPERP1 * std::sin(ALPHA1);
1285 G4double VY1_IMF = VPERP1 * std::cos(ALPHA1);
1289
1290 G4double EEIMFP = EEBU * AFBU / (AFBU + AIMFBU);
1291 G4double EEIMF = EEBU * AIMFBU / (AFBU + AIMFBU);
1292
1293
1295 0.40 * 931.490 * 1.160 * 1.160 * (std::pow(AIMFBU, 5.0 / 3.0) + std::pow(AFBU, 5.0 / 3.0)) +
1296 931.490 * 1.160 * 1.160 * AIMFBU * AFBU / (AIMFBU + AFBU) *
1297 (std::pow(AIMFBU, 1. / 3.) + std::pow(AFBU, 1. / 3.)) *
1298 (std::pow(AIMFBU, 1. / 3.) + std::pow(AFBU, 1. / 3.));
1299
1301 BU_TAB[i][9] * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(AFBU, 5.0 / 3.0) / IINERTTOT;
1303 BU_TAB[i][9] * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(AIMFBU, 5.0 / 3.0) / IINERTTOT;
1304
1305
1306
1307
1308
1309
1311 VX1_IMF, VY1_IMF, VZ1_IMF, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VXOUT, &VYOUT, &VZOUT);
1312 BU_TAB[i][4] = VXOUT;
1313 BU_TAB[i][5] = VYOUT;
1314 BU_TAB[i][6] = VZOUT;
1315
1316 G4double vx1ev_imf = 0., vy1ev_imf = 0., vz1ev_imf = 0., zdummy = 0., adummy = 0., tkedummy = 0.,
1317 jprf1 = 0.;
1318
1319
1322 G4double pbH = (AFBU - ZFBU) / (AFBU - ZFBU + AIMFBU - ZIMFBU);
1323 for (
G4int j = 0; j < nbl; j++)
1324 {
1326 {
1327 NbLamH++;
1328 }
1329 else
1330 {
1331 NbLamimf++;
1332 }
1333 }
1334
1336 AFBU,
1337 &EEIMFP,
1338 JPRFHEAVY,
1339 &ZFFBU,
1340 &AFFBU,
1341 &mtota,
1342 &vz1ev_imf,
1343 &vx1ev_imf,
1344 &vy1ev_imf,
1345 &FFBU1,
1346 &FIMFBU1,
1347 &zdummy,
1348 &adummy,
1349 &tkedummy,
1350 &jprf1,
1351 &inttype,
1352 &inum,
1353 EV_TEMP,
1354 &IEV_TAB_TEMP,
1355 &NbLamH);
1356
1357 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++)
1358 {
1359 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1360 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1361 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1362
1363
1364
1365
1366
1368 BU_TAB[i][5],
1369 BU_TAB[i][6],
1370 EV_TEMP[IJ][2],
1371 EV_TEMP[IJ][3],
1372 EV_TEMP[IJ][4],
1373 &VXOUT,
1374 &VYOUT,
1375 &VZOUT);
1376 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1377 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1378 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1379 }
1380 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1381
1382 BU_TAB[i][7] =
dint(ZFFBU);
1383 BU_TAB[i][8] =
dint(AFFBU);
1384 BU_TAB[i][11] = NbLamH;
1385
1386
1387
1388
1390 vy1ev_imf,
1391 vz1ev_imf,
1392 BU_TAB[i][4],
1393 BU_TAB[i][5],
1394 BU_TAB[i][6],
1395 &VXOUT,
1396 &VYOUT,
1397 &VZOUT);
1398 BU_TAB[i][4] = VXOUT;
1399 BU_TAB[i][5] = VYOUT;
1400 BU_TAB[i][6] = VZOUT;
1401
1404 opt->optimfallowed = 0;
1405 fiss->ifis = 0;
1406
1407 G4double zffimf, affimf, zdummy1, adummy1, tkedummy1, jprf2, vx2ev_imf, vy2ev_imf, vz2ev_imf;
1408
1410 AIMFBU,
1411 &EEIMF,
1412 JPRFLIGHT,
1413 &zffimf,
1414 &affimf,
1415 &mtota,
1416 &vz2ev_imf,
1417 &vx2ev_imf,
1418 &vy2ev_imf,
1419 &FFBU2,
1420 &FIMFBU2,
1421 &zdummy1,
1422 &adummy1,
1423 &tkedummy1,
1424 &jprf2,
1425 &inttype,
1426 &inum,
1427 EV_TEMP,
1428 &IEV_TAB_TEMP,
1429 &NbLamimf);
1430
1431 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++)
1432 {
1433 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1434 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1435 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1436
1437
1438
1439
1440
1441
1443 BU_TAB[i][5],
1444 BU_TAB[i][6],
1445 EV_TEMP[IJ][2],
1446 EV_TEMP[IJ][3],
1447 EV_TEMP[IJ][4],
1448 &VXOUT,
1449 &VYOUT,
1450 &VZOUT);
1451 lorentz_boost(VX2_IMF, VY2_IMF, VZ2_IMF, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1452 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
1453 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
1454 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
1455 }
1456 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1457
1458 BU_TAB[IMULTBU + ILOOPBU][0] = BU_TAB[i][0];
1459 BU_TAB[IMULTBU + ILOOPBU][1] = BU_TAB[i][1];
1460 BU_TAB[IMULTBU + ILOOPBU][2] = BU_TAB[i][2];
1461 BU_TAB[IMULTBU + ILOOPBU][3] = BU_TAB[i][3];
1462 BU_TAB[IMULTBU + ILOOPBU][7] =
dint(zffimf);
1463 BU_TAB[IMULTBU + ILOOPBU][8] =
dint(affimf);
1464 BU_TAB[IMULTBU + ILOOPBU][11] = NbLamimf;
1465
1467 VX2_IMF, VY2_IMF, VZ2_IMF, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VXOUT, &VYOUT, &VZOUT);
1468 lorentz_boost(vx2ev_imf, vy2ev_imf, vz2ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1469 BU_TAB[IMULTBU + ILOOPBU][4] = VX2OUT;
1470 BU_TAB[IMULTBU + ILOOPBU][5] = VY2OUT;
1471 BU_TAB[IMULTBU + ILOOPBU][6] = VZ2OUT;
1472 ILOOPBU = ILOOPBU + 1;
1473 }
1474 }
1475 else
1476 {
1477
1478
1479
1480
1481 BU_TAB[i][7] = BU_TAB[i][0];
1482 BU_TAB[i][8] = BU_TAB[i][1];
1483
1484
1485
1486 BU_TAB[i][11] = Nblamb[i];
1487 }
1488 }
1489
1490 IMULTBU = IMULTBU + ILOOPBU;
1491
1492
1493
1494 INEWLOOP = 0;
1495 ABU_SUM = 0.0;
1496 ZBU_SUM = 0.0;
1497
1498 for (
G4int i = 0; i < IMULTBU; i++)
1499 {
1500 ABU_SUM = ABU_SUM + BU_TAB[i][8];
1501 ZBU_SUM = ZBU_SUM + BU_TAB[i][7];
1504 &afpnew,
1505 &zfpnew,
1506 IOUNSTABLE,
1507 BU_TAB[i][4],
1508 BU_TAB[i][5],
1509 BU_TAB[i][6],
1510 &VP1X,
1511 &VP1Y,
1512 &VP1Z,
1513 BU_TAB_TEMP1,
1514 &ILOOP);
1515
1516
1517
1518
1519
1520
1521 if (IOUNSTABLE > 0)
1522 {
1523
1524 ABU_SUM = ABU_SUM +
G4double(afpnew) - BU_TAB[i][8];
1525 ZBU_SUM = ZBU_SUM +
G4double(zfpnew) - BU_TAB[i][7];
1528 BU_TAB[i][4] = VP1X;
1529 BU_TAB[i][5] = VP1Y;
1530 BU_TAB[i][6] = VP1Z;
1531
1532
1533 for (
G4int IJ = 0; IJ < ILOOP; IJ++)
1534 {
1535 BU_TAB[IMULTBU + INEWLOOP + IJ][7] = BU_TAB_TEMP1[IJ][0];
1536 BU_TAB[IMULTBU + INEWLOOP + IJ][8] = BU_TAB_TEMP1[IJ][1];
1537 BU_TAB[IMULTBU + INEWLOOP + IJ][4] = BU_TAB_TEMP1[IJ][2];
1538 BU_TAB[IMULTBU + INEWLOOP + IJ][5] = BU_TAB_TEMP1[IJ][3];
1539 BU_TAB[IMULTBU + INEWLOOP + IJ][6] = BU_TAB_TEMP1[IJ][4];
1540 BU_TAB[IMULTBU + INEWLOOP + IJ][2] = 0.0;
1541 BU_TAB[IMULTBU + INEWLOOP + IJ][3] = 0.0;
1542 BU_TAB[IMULTBU + INEWLOOP + IJ][0] = BU_TAB[i][0];
1543 BU_TAB[IMULTBU + INEWLOOP + IJ][1] = BU_TAB[i][1];
1544 BU_TAB[IMULTBU + INEWLOOP + IJ][11] = BU_TAB[i][11];
1545 ABU_SUM = ABU_SUM + BU_TAB[IMULTBU + INEWLOOP + IJ][8];
1546 ZBU_SUM = ZBU_SUM + BU_TAB[IMULTBU + INEWLOOP + IJ][7];
1547 }
1548
1549 INEWLOOP = INEWLOOP + ILOOP;
1550 }
1551 }
1552
1553
1554 IMULTBU = IMULTBU + INEWLOOP;
1555
1556
1557 lorentz_boost(VX_incl, VY_incl, VZ_incl, VX_PREF, VY_PREF, VZ_PREF, &VXOUT, &VYOUT, &VZOUT);
1558 VX_PREF = VXOUT;
1559 VY_PREF = VYOUT;
1560 VZ_PREF = VZOUT;
1561
1562 for (
G4int i = 0; i < IMULTBU; i++)
1563 {
1564 lorentz_boost(VX_incl, VY_incl, VZ_incl, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VXOUT, &VYOUT, &VZOUT);
1565 BU_TAB[i][4] = VXOUT;
1566 BU_TAB[i][5] = VYOUT;
1567 BU_TAB[i][6] = VZOUT;
1568 }
1569 for (
G4int i = 0; i < IEV_TAB; i++)
1570 {
1571 lorentz_boost(VX_incl, VY_incl, VZ_incl, EV_TAB[i][2], EV_TAB[i][3], EV_TAB[i][4], &VXOUT, &VYOUT, &VZOUT);
1572 EV_TAB[i][2] = VXOUT;
1573 EV_TAB[i][3] = VYOUT;
1574 EV_TAB[i][4] = VZOUT;
1575 }
1576 if (IMULTBU > 200)
1577 std::cout << "IMULTBU>200 " << IMULTBU << std::endl;
1578 delete[] problamb;
1579 delete[] Nblamb;
1580 }
1581
1582mult7777:
1583
1584
1587
1588 if (IMULTIFR == 0)
1589 {
1590
1591
1592 VX_PREF = VX_incl;
1593 VY_PREF = VY_incl;
1594 VZ_PREF = VZ_incl;
1595 }
1596
1597 if (IMULTIFR == 1)
1598 {
1599 NbLam0 = NbLamprf;
1600 }
1601
1602
1603
1604 opt->optimfallowed = 1;
1605 fiss->ifis = 1;
1606 fimf = 0;
1607 ff = 0;
1608
1609
1610
1611 if (zprfp <= 2 && zprfp < aprfp)
1612 {
1613 zf = zprf;
1614 af = aprf;
1615 ee = 0.0;
1616 ff = 0;
1617 fimf = 0;
1618 ftype = 0;
1619 aimf = 0.0;
1620 zimf = 0.0;
1621 tkeimf = 0.0;
1622 vx_eva = 0.0;
1623 vy_eva = 0.0;
1624 vz_eva = 0.0;
1625 jprf0 = jprf;
1626 goto a1972;
1627 }
1628
1629
1630 if (zprfp <= 2 && zprfp == aprfp)
1631 {
1633 zprfp,
1634 &afpnew,
1635 &zfpnew,
1636 IOUNSTABLE,
1637 VX_PREF,
1638 VY_PREF,
1639 VZ_PREF,
1640 &VP1X,
1641 &VP1Y,
1642 &VP1Z,
1643 EV_TAB_TEMP,
1644 &ILOOP);
1647 VX_PREF = VP1X;
1648 VY_PREF = VP1Y;
1649 VZ_PREF = VP1Z;
1650 for (
G4int I = 0; I < ILOOP; I++)
1651 {
1652 for (
G4int IJ = 0; IJ < 6; IJ++)
1653 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1654 }
1655 IEV_TAB = IEV_TAB + ILOOP;
1656 ee = 0.0;
1657 ff = 0;
1658 fimf = 0;
1659 ftype = 0;
1660 aimf = 0.0;
1661 zimf = 0.0;
1662 tkeimf = 0.0;
1663 vx_eva = 0.0;
1664 vy_eva = 0.0;
1665 vz_eva = 0.0;
1666 jprf0 = jprf;
1667 goto a1972;
1668 }
1669
1670
1671 if (zprfp == aprfp)
1672 {
1674 zprfp,
1675 &afpnew,
1676 &zfpnew,
1677 IOUNSTABLE,
1678 VX_PREF,
1679 VY_PREF,
1680 VZ_PREF,
1681 &VP1X,
1682 &VP1Y,
1683 &VP1Z,
1684 EV_TAB_TEMP,
1685 &ILOOP);
1688 VX_PREF = VP1X;
1689 VY_PREF = VP1Y;
1690 VZ_PREF = VP1Z;
1691 for (
G4int I = 0; I < ILOOP; I++)
1692 {
1693 for (
G4int IJ = 0; IJ < 6; IJ++)
1694 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1695 }
1696 IEV_TAB = IEV_TAB + ILOOP;
1697 ee = 0.0;
1698 ff = 0;
1699 fimf = 0;
1700 ftype = 0;
1701 aimf = 0.0;
1702 zimf = 0.0;
1703 tkeimf = 0.0;
1704 vx_eva = 0.0;
1705 vy_eva = 0.0;
1706 vz_eva = 0.0;
1707 jprf0 = jprf;
1708 goto a1972;
1709 }
1710
1712 aprf,
1713 &ee,
1714 jprf,
1715 &zf,
1716 &af,
1717 &mtota,
1718 &vz_eva,
1719 &vx_eva,
1720 &vy_eva,
1721 &ff,
1722 &fimf,
1723 &zimf,
1724 &aimf,
1725 &tkeimf,
1726 &jprf0,
1727 &inttype,
1728 &inum,
1729 EV_TEMP,
1730 &IEV_TAB_TEMP,
1731 &NbLam0);
1732
1733 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++)
1734 {
1735 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1736 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1737 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1738
1739
1740
1741
1742
1744 VX_PREF, VY_PREF, VZ_PREF, EV_TEMP[IJ][2], EV_TEMP[IJ][3], EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1745 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1746 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1747 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1748 }
1749 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1750
1751a1972:
1752
1753
1754 lorentz_boost(VX_PREF, VY_PREF, VZ_PREF, vx_eva, vy_eva, vz_eva, &VXOUT, &VYOUT, &VZOUT);
1755 V_CM[0] = VXOUT;
1756 V_CM[1] = VYOUT;
1757 V_CM[2] = VZOUT;
1758
1759 if (ff == 0 && fimf == 0)
1760 {
1761
1762 ftype = 0;
1765 SFP1 = NbLam0;
1766 AFPIMF = 0;
1767 ZFPIMF = 0;
1768 SFPIMF = 0;
1769 ZFP2 = 0;
1770 AFP2 = 0;
1771 SFP2 = 0;
1772 VFP1_CM[0] = V_CM[0];
1773 VFP1_CM[1] = V_CM[1];
1774 VFP1_CM[2] = V_CM[2];
1775 for (
G4int j = 0; j < 3; j++)
1776 {
1777 VIMF_CM[j] = 0.0;
1778 VFP2_CM[j] = 0.0;
1779 }
1780 }
1781
1782 if (ff == 1 && fimf == 0)
1783 ftype = 1;
1784 if (ff == 0 && fimf == 1)
1785 ftype = 2;
1786
1787
1788
1789
1790
1791
1792
1793 if (ftype == 1)
1794 {
1795 varntp->kfis = 1;
1796 if (NbLam0 > 0)
1797 varntp->kfis = 20;
1798
1799
1800 G4int IEV_TAB_FIS = 0, imode = 0;
1801
1802 G4double vx1_fission = 0., vy1_fission = 0., vz1_fission = 0.;
1803 G4double vx2_fission = 0., vy2_fission = 0., vz2_fission = 0.;
1804 G4double vx_eva_sc = 0., vy_eva_sc = 0., vz_eva_sc = 0.;
1805
1807 zf,
1808 ee,
1809 jprf0,
1810 &vx1_fission,
1811 &vy1_fission,
1812 &vz1_fission,
1813 &vx2_fission,
1814 &vy2_fission,
1815 &vz2_fission,
1816 &ZFP1,
1817 &AFP1,
1818 &SFP1,
1819 &ZFP2,
1820 &AFP2,
1821 &SFP2,
1822 &imode,
1823 &vx_eva_sc,
1824 &vy_eva_sc,
1825 &vz_eva_sc,
1826 EV_TEMP,
1827 &IEV_TAB_FIS,
1828 &NbLam0);
1829
1830 for (
G4int IJ = 0; IJ < IEV_TAB_FIS; IJ++)
1831 {
1832 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1833 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1834 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1835
1836
1837
1838
1839
1841 V_CM[0], V_CM[1], V_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3], EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1842 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1843 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1844 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1845 }
1846 IEV_TAB = IEV_TAB + IEV_TAB_FIS;
1847
1848
1849
1850
1851 AFPIMF = 0;
1852 ZFPIMF = 0;
1853 SFPIMF = 0;
1854
1855
1856
1857
1858
1859
1860
1861
1862 lorentz_boost(vx1_fission, vy1_fission, vz1_fission, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
1863 lorentz_boost(vx_eva_sc, vy_eva_sc, vz_eva_sc, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1864 VFP1_CM[0] = VX2OUT;
1865 VFP1_CM[1] = VY2OUT;
1866 VFP1_CM[2] = VZ2OUT;
1867
1868
1869
1870
1871
1872
1873 lorentz_boost(vx2_fission, vy2_fission, vz2_fission, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
1874 lorentz_boost(vx_eva_sc, vy_eva_sc, vz_eva_sc, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1875 VFP2_CM[0] = VX2OUT;
1876 VFP2_CM[1] = VY2OUT;
1877 VFP2_CM[2] = VZ2OUT;
1878
1879
1880
1881
1882 }
1883 else if (ftype == 2)
1884 {
1885
1886
1889 opt->optimfallowed = 1;
1890 fiss->ifis = 1;
1891
1894 G4double pbH = (af - zf) / (af - zf + aimf - zimf);
1895
1896 for (
G4int i = 0; i < NbLam0; i++)
1897 {
1899 {
1900 NbLamH++;
1901 }
1902 else
1903 {
1904 NbLamimf++;
1905 }
1906 }
1907
1908
1909 G4double EkinR1 = tkeimf * aimf / (af + aimf);
1910 G4double EkinR2 = tkeimf * af / (af + aimf);
1911 G4double V1 = std::sqrt(EkinR1 / af) * 1.3887;
1912 G4double V2 = std::sqrt(EkinR2 / aimf) * 1.3887;
1914 G4double VPERP1 = std::sqrt(
V1 *
V1 - VZ1_IMF * VZ1_IMF);
1916 G4double VX1_IMF = VPERP1 * std::sin(ALPHA1);
1917 G4double VY1_IMF = VPERP1 * std::cos(ALPHA1);
1921
1922 G4double EEIMFP = ee * af / (af + aimf);
1923 G4double EEIMF = ee * aimf / (af + aimf);
1924
1925
1926 G4double IINERTTOT = 0.40 * 931.490 * 1.160 * 1.160 * (std::pow(aimf, 5.0 / 3.0) + std::pow(af, 5.0 / 3.0)) +
1927 931.490 * 1.160 * 1.160 * aimf * af / (aimf + af) *
1928 (std::pow(aimf, 1. / 3.) + std::pow(af, 1. / 3.)) *
1929 (std::pow(aimf, 1. / 3.) + std::pow(af, 1. / 3.));
1930
1931 G4double JPRFHEAVY = jprf0 * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(af, 5.0 / 3.0) / IINERTTOT;
1932 G4double JPRFLIGHT = jprf0 * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(aimf, 5.0 / 3.0) / IINERTTOT;
1933 if (af < 2.0)
1934 std::cout << "RN117-4,AF,ZF,EE,JPRFheavy" << std::endl;
1935
1936 G4double vx1ev_imf = 0., vy1ev_imf = 0., vz1ev_imf = 0., zdummy = 0., adummy = 0., tkedummy = 0., jprf1 = 0.;
1937
1939 af,
1940 &EEIMFP,
1941 JPRFHEAVY,
1942 &zff,
1943 &aff,
1944 &mtota,
1945 &vz1ev_imf,
1946 &vx1ev_imf,
1947 &vy1ev_imf,
1948 &FF11,
1949 &FIMF11,
1950 &zdummy,
1951 &adummy,
1952 &tkedummy,
1953 &jprf1,
1954 &inttype,
1955 &inum,
1956 EV_TEMP,
1957 &IEV_TAB_TEMP,
1958 &NbLamH);
1959
1960 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++)
1961 {
1962 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1963 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1964 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1965
1966
1967
1968
1969
1971 V_CM[0], V_CM[1], V_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3], EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1972 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1973 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
1974 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
1975 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
1976 }
1977 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1978
1979
1982 opt->optimfallowed = 0;
1983 fiss->ifis = 0;
1984
1985
1986 G4double zffimf, affimf, zdummy1 = 0., adummy1 = 0., tkedummy1 = 0., jprf2, vx2ev_imf, vy2ev_imf, vz2ev_imf;
1987
1989 aimf,
1990 &EEIMF,
1991 JPRFLIGHT,
1992 &zffimf,
1993 &affimf,
1994 &mtota,
1995 &vz2ev_imf,
1996 &vx2ev_imf,
1997 &vy2ev_imf,
1998 &FF22,
1999 &FIMF22,
2000 &zdummy1,
2001 &adummy1,
2002 &tkedummy1,
2003 &jprf2,
2004 &inttype,
2005 &inum,
2006 EV_TEMP,
2007 &IEV_TAB_TEMP,
2008 &NbLamimf);
2009
2010 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++)
2011 {
2012 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
2013 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
2014 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
2015
2016
2017
2018
2019
2021 V_CM[0], V_CM[1], V_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3], EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
2022 lorentz_boost(VX2_IMF, VY2_IMF, VZ2_IMF, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2023 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
2024 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
2025 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
2026 }
2027 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
2028
2029
2032 SFPIMF = NbLamimf;
2033
2034
2035
2036
2037
2038
2039
2040 lorentz_boost(VX2_IMF, VY2_IMF, VZ2_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
2041 lorentz_boost(vx2ev_imf, vy2ev_imf, vz2ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2042 VIMF_CM[0] = VX2OUT;
2043 VIMF_CM[1] = VY2OUT;
2044 VIMF_CM[2] = VZ2OUT;
2045
2046
2047
2048
2049 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
2050 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2051 VFP1_CM[0] = VX2OUT;
2052 VFP1_CM[1] = VY2OUT;
2053 VFP1_CM[2] = VZ2OUT;
2054
2055 if (FF11 == 0 && FIMF11 == 0)
2056 {
2057
2060 SFP1 = NbLamH;
2061 ZFP2 = 0;
2062 AFP2 = 0;
2063 SFP2 = 0;
2064 ftype = 2;
2067 SFPIMF = NbLamimf;
2068 for (
G4int I = 0; I < 3; I++)
2069 VFP2_CM[I] = 0.0;
2070 }
2071 else if (FF11 == 1 && FIMF11 == 0)
2072 {
2073
2074 varntp->kfis = 1;
2075 if (NbLam0 > 0)
2076 varntp->kfis = 20;
2077
2078 opt->optimfallowed = 0;
2079 fiss->ifis = 0;
2080
2081 zf = zff;
2082 af = aff;
2083 ee = EEIMFP;
2084
2085 ftype = 21;
2086
2087 G4int IEV_TAB_FIS = 0, imode = 0;
2088
2089 G4double vx1_fission = 0., vy1_fission = 0., vz1_fission = 0.;
2090 G4double vx2_fission = 0., vy2_fission = 0., vz2_fission = 0.;
2091 G4double vx_eva_sc = 0., vy_eva_sc = 0., vz_eva_sc = 0.;
2092
2094 zf,
2095 ee,
2096 jprf1,
2097 &vx1_fission,
2098 &vy1_fission,
2099 &vz1_fission,
2100 &vx2_fission,
2101 &vy2_fission,
2102 &vz2_fission,
2103 &ZFP1,
2104 &AFP1,
2105 &SFP1,
2106 &ZFP2,
2107 &AFP2,
2108 &SFP2,
2109 &imode,
2110 &vx_eva_sc,
2111 &vy_eva_sc,
2112 &vz_eva_sc,
2113 EV_TEMP,
2114 &IEV_TAB_FIS,
2115 &NbLamH);
2116
2117 for (int IJ = 0; IJ < IEV_TAB_FIS; IJ++)
2118 {
2119 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
2120 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
2121 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
2122
2123
2124
2125
2126
2128 VFP1_CM[1],
2129 VFP1_CM[2],
2130 EV_TEMP[IJ][2],
2131 EV_TEMP[IJ][3],
2132 EV_TEMP[IJ][4],
2133 &VXOUT,
2134 &VYOUT,
2135 &VZOUT);
2136 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
2137 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
2138 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
2139 }
2140 IEV_TAB = IEV_TAB + IEV_TAB_FIS;
2141
2142
2143
2144
2145
2146
2147
2148
2149
2150
2151
2152 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
2153 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2154 lorentz_boost(vx1_fission, vy1_fission, vz1_fission, VX2OUT, VY2OUT, VZ2OUT, &VXOUT, &VYOUT, &VZOUT);
2155 lorentz_boost(vx_eva_sc, vy_eva_sc, vz_eva_sc, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2156 VFP1_CM[0] = VX2OUT;
2157 VFP1_CM[1] = VY2OUT;
2158 VFP1_CM[2] = VZ2OUT;
2159
2160
2161
2162
2163
2164
2165
2166
2167 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
2168 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2169 lorentz_boost(vx2_fission, vy2_fission, vz2_fission, VX2OUT, VY2OUT, VZ2OUT, &VXOUT, &VYOUT, &VZOUT);
2170 lorentz_boost(vx_eva_sc, vy_eva_sc, vz_eva_sc, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2171 VFP2_CM[0] = VX2OUT;
2172 VFP2_CM[1] = VY2OUT;
2173 VFP2_CM[2] = VZ2OUT;
2174 }
2175 else if (FF11 == 0 && FIMF11 == 1)
2176 {
2177
2178
2179 opt->optimfallowed = 0;
2180 fiss->ifis = 0;
2181
2182 zf = zff;
2183 af = aff;
2184 ee = EEIMFP;
2185 aimf = adummy;
2186 zimf = zdummy;
2187 tkeimf = tkedummy;
2188 FF11 = 0;
2189 FIMF11 = 0;
2190 ftype = 22;
2191
2193 G4int NbLamimf1 = 0;
2194 G4double pbH1 = (af - zf) / (af - zf + aimf - zimf);
2195 for (
G4int i = 0; i < NbLamH; i++)
2196 {
2198 {
2199 NbLamH1++;
2200 }
2201 else
2202 {
2203 NbLamimf1++;
2204 }
2205 }
2206
2207
2208 EkinR1 = tkeimf * aimf / (af + aimf);
2209 EkinR2 = tkeimf * af / (af + aimf);
2210 V1 = std::sqrt(EkinR1 / af) * 1.3887;
2211 V2 = std::sqrt(EkinR2 / aimf) * 1.3887;
2213 VPERP1 = std::sqrt(
V1 *
V1 - VZ1_IMFS * VZ1_IMFS);
2215 G4double VX1_IMFS = VPERP1 * std::sin(ALPHA1);
2216 G4double VY1_IMFS = VPERP1 * std::cos(ALPHA1);
2220
2221 EEIMFP = ee * af / (af + aimf);
2222 EEIMF = ee * aimf / (af + aimf);
2223
2224
2225 IINERTTOT = 0.40 * 931.490 * 1.160 * 1.160 * (std::pow(aimf, 5.0 / 3.0) + std::pow(af, 5.0 / 3.0)) +
2226 931.490 * 1.160 * 1.160 * aimf * af / (aimf + af) *
2227 (std::pow(aimf, 1. / 3.) + std::pow(af, 1. / 3.)) *
2228 (std::pow(aimf, 1. / 3.) + std::pow(af, 1. / 3.));
2229
2230 JPRFHEAVY = jprf1 * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(af, 5.0 / 3.0) / IINERTTOT;
2231 JPRFLIGHT = jprf1 * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(aimf, 5.0 / 3.0) / IINERTTOT;
2232
2233 G4double zffs = 0., affs = 0., vx1ev_imfs = 0., vy1ev_imfs = 0., vz1ev_imfs = 0., jprf3 = 0.;
2234
2236 af,
2237 &EEIMFP,
2238 JPRFHEAVY,
2239 &zffs,
2240 &affs,
2241 &mtota,
2242 &vz1ev_imfs,
2243 &vx1ev_imfs,
2244 &vy1ev_imfs,
2245 &FF11,
2246 &FIMF11,
2247 &zdummy,
2248 &adummy,
2249 &tkedummy,
2250 &jprf3,
2251 &inttype,
2252 &inum,
2253 EV_TEMP,
2254 &IEV_TAB_TEMP,
2255 &NbLamH1);
2256
2257 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++)
2258 {
2259 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
2260 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
2261 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
2262
2263
2264
2265
2266
2268 VFP1_CM[1],
2269 VFP1_CM[2],
2270 EV_TEMP[IJ][2],
2271 EV_TEMP[IJ][3],
2272 EV_TEMP[IJ][4],
2273 &VXOUT,
2274 &VYOUT,
2275 &VZOUT);
2276 lorentz_boost(vx1ev_imfs, vy1ev_imfs, vz1ev_imfs, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2277 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
2278 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
2279 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
2280 }
2281 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
2282
2283
2284 opt->optimfallowed = 0;
2285 fiss->ifis = 0;
2286
2287 FF22 = 0;
2288 FIMF22 = 0;
2289
2290 G4double zffimfs = 0., affimfs = 0., vx2ev_imfs = 0., vy2ev_imfs = 0., vz2ev_imfs = 0., jprf4 = 0.;
2291
2293 aimf,
2294 &EEIMF,
2295 JPRFLIGHT,
2296 &zffimfs,
2297 &affimfs,
2298 &mtota,
2299 &vz2ev_imfs,
2300 &vx2ev_imfs,
2301 &vy2ev_imfs,
2302 &FF22,
2303 &FIMF22,
2304 &zdummy1,
2305 &adummy1,
2306 &tkedummy1,
2307 &jprf4,
2308 &inttype,
2309 &inum,
2310 EV_TEMP,
2311 &IEV_TAB_TEMP,
2312 &NbLamimf1);
2313
2314 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++)
2315 {
2316 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
2317 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
2318 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
2319
2320
2321
2322
2323
2325 VFP1_CM[1],
2326 VFP1_CM[2],
2327 EV_TEMP[IJ][2],
2328 EV_TEMP[IJ][3],
2329 EV_TEMP[IJ][4],
2330 &VXOUT,
2331 &VYOUT,
2332 &VZOUT);
2333 lorentz_boost(vx2ev_imfs, vy2ev_imfs, vz2ev_imfs, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2334 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
2335 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
2336 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
2337 }
2338 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
2339
2342 SFP1 = NbLamH1;
2345 SFP2 = NbLamimf1;
2346
2347
2348
2349
2350
2351
2352 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
2353 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2354 lorentz_boost(VX1_IMFS, VY1_IMFS, VZ1_IMFS, VX2OUT, VY2OUT, VZ2OUT, &VXOUT, &VYOUT, &VZOUT);
2355 lorentz_boost(vx1ev_imfs, vy1ev_imfs, vz1ev_imfs, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2356 VFP1_CM[0] = VX2OUT;
2357 VFP1_CM[1] = VY2OUT;
2358 VFP1_CM[2] = VZ2OUT;
2359
2360
2361
2362
2363
2364
2365 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
2366 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2367 lorentz_boost(VX2_IMFS, VY2_IMFS, VZ2_IMFS, VX2OUT, VY2OUT, VZ2OUT, &VXOUT, &VYOUT, &VZOUT);
2368 lorentz_boost(vx2ev_imfs, vy2ev_imfs, vz2ev_imfs, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
2369 VFP2_CM[0] = VX2OUT;
2370 VFP2_CM[1] = VY2OUT;
2371 VFP2_CM[2] = VZ2OUT;
2372 }
2373 }
2374
2375
2376 if (ftype != 1 && ftype != 21)
2377 {
2378
2379
2380 IOUNSTABLE = 0;
2381
2383 ZFP1,
2384 &afpnew,
2385 &zfpnew,
2386 IOUNSTABLE,
2387 VFP1_CM[0],
2388 VFP1_CM[1],
2389 VFP1_CM[2],
2390 &VP1X,
2391 &VP1Y,
2392 &VP1Z,
2393 EV_TAB_TEMP,
2394 &ILOOP);
2395
2396 if (IOUNSTABLE == 1)
2397 {
2398 AFP1 = afpnew;
2399 ZFP1 = zfpnew;
2400 VFP1_CM[0] = VP1X;
2401 VFP1_CM[1] = VP1Y;
2402 VFP1_CM[2] = VP1Z;
2403 for (
G4int I = 0; I < ILOOP; I++)
2404 {
2405 for (
G4int IJ = 0; IJ < 5; IJ++)
2406 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2407 }
2408 IEV_TAB = IEV_TAB + ILOOP;
2409 }
2410
2411 if (ftype > 1)
2412 {
2413 IOUNSTABLE = 0;
2414
2416 ZFPIMF,
2417 &afpnew,
2418 &zfpnew,
2419 IOUNSTABLE,
2420 VIMF_CM[0],
2421 VIMF_CM[1],
2422 VIMF_CM[2],
2423 &VP1X,
2424 &VP1Y,
2425 &VP1Z,
2426 EV_TAB_TEMP,
2427 &ILOOP);
2428
2429 if (IOUNSTABLE == 1)
2430 {
2431 AFPIMF = afpnew;
2432 ZFPIMF = zfpnew;
2433 VIMF_CM[0] = VP1X;
2434 VIMF_CM[1] = VP1Y;
2435 VIMF_CM[2] = VP1Z;
2436 for (
G4int I = 0; I < ILOOP; I++)
2437 {
2438 for (
G4int IJ = 0; IJ < 5; IJ++)
2439 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2440 }
2441 IEV_TAB = IEV_TAB + ILOOP;
2442 }
2443
2444 if (ftype > 2)
2445 {
2446 IOUNSTABLE = 0;
2447
2449 ZFP2,
2450 &afpnew,
2451 &zfpnew,
2452 IOUNSTABLE,
2453 VFP2_CM[0],
2454 VFP2_CM[1],
2455 VFP2_CM[2],
2456 &VP1X,
2457 &VP1Y,
2458 &VP1Z,
2459 EV_TAB_TEMP,
2460 &ILOOP);
2461
2462 if (IOUNSTABLE == 1)
2463 {
2464 AFP2 = afpnew;
2465 ZFP2 = zfpnew;
2466 VFP2_CM[0] = VP1X;
2467 VFP2_CM[1] = VP1Y;
2468 VFP2_CM[2] = VP1Z;
2469 for (
G4int I = 0; I < ILOOP; I++)
2470 {
2471 for (
G4int IJ = 0; IJ < 5; IJ++)
2472 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2473 }
2474 IEV_TAB = IEV_TAB + ILOOP;
2475 }
2476 }
2477 }
2478 }
2479
2480
2481 if (ftype == 1 || ftype == 21)
2482 {
2483
2484 IOUNSTABLE = 0;
2485
2487 ZFP1,
2488 &afpnew,
2489 &zfpnew,
2490 IOUNSTABLE,
2491 VFP1_CM[0],
2492 VFP1_CM[1],
2493 VFP1_CM[2],
2494 &VP1X,
2495 &VP1Y,
2496 &VP1Z,
2497 EV_TAB_TEMP,
2498 &ILOOP);
2499
2500 if (IOUNSTABLE == 1)
2501 {
2502 AFP1 = afpnew;
2503 ZFP1 = zfpnew;
2504 VFP1_CM[0] = VP1X;
2505 VFP1_CM[1] = VP1Y;
2506 VFP1_CM[2] = VP1Z;
2507 for (
G4int I = 0; I < ILOOP; I++)
2508 {
2509 for (
G4int IJ = 0; IJ < 5; IJ++)
2510 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2511 }
2512 IEV_TAB = IEV_TAB + ILOOP;
2513 }
2514
2515 IOUNSTABLE = 0;
2516
2518 ZFP2,
2519 &afpnew,
2520 &zfpnew,
2521 IOUNSTABLE,
2522 VFP2_CM[0],
2523 VFP2_CM[1],
2524 VFP2_CM[2],
2525 &VP1X,
2526 &VP1Y,
2527 &VP1Z,
2528 EV_TAB_TEMP,
2529 &ILOOP);
2530
2531 if (IOUNSTABLE == 1)
2532 {
2533 AFP2 = afpnew;
2534 ZFP2 = zfpnew;
2535 VFP2_CM[0] = VP1X;
2536 VFP2_CM[1] = VP1Y;
2537 VFP2_CM[2] = VP1Z;
2538 for (
G4int I = 0; I < ILOOP; I++)
2539 {
2540 for (
G4int IJ = 0; IJ < 5; IJ++)
2541 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2542 }
2543 IEV_TAB = IEV_TAB + ILOOP;
2544 }
2545
2546 if (ftype == 21)
2547 {
2548 IOUNSTABLE = 0;
2549
2551 ZFPIMF,
2552 &afpnew,
2553 &zfpnew,
2554 IOUNSTABLE,
2555 VIMF_CM[0],
2556 VIMF_CM[1],
2557 VIMF_CM[2],
2558 &VP1X,
2559 &VP1Y,
2560 &VP1Z,
2561 EV_TAB_TEMP,
2562 &ILOOP);
2563
2564 if (IOUNSTABLE == 1)
2565 {
2566 AFPIMF = afpnew;
2567 ZFPIMF = zfpnew;
2568 VIMF_CM[0] = VP1X;
2569 VIMF_CM[1] = VP1Y;
2570 VIMF_CM[2] = VP1Z;
2571 for (
G4int I = 0; I < ILOOP; I++)
2572 {
2573 for (
G4int IJ = 0; IJ < 5; IJ++)
2574 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2575 }
2576 IEV_TAB = IEV_TAB + ILOOP;
2577 }
2578 }
2579 }
2580
2581
2582 if ((ftype == 1 || ftype == 21) && (AFP2 <= 0 || AFP1 <= 0 || ZFP2 <= 0 || ZFP1 <= 0))
2583 {
2584 std::cout << "ZFP1:" << ZFP1 << std::endl;
2585 std::cout << "AFP1:" << AFP1 << std::endl;
2586 std::cout << "ZFP2:" << ZFP2 << std::endl;
2587 std::cout << "AFP2:" << AFP2 << std::endl;
2588 }
2589
2590
2591 EV_TAB[IEV_TAB][0] = ZFP1;
2592 EV_TAB[IEV_TAB][1] = AFP1;
2593 EV_TAB[IEV_TAB][5] = SFP1;
2594 EV_TAB[IEV_TAB][2] = VFP1_CM[0];
2595 EV_TAB[IEV_TAB][3] = VFP1_CM[1];
2596 EV_TAB[IEV_TAB][4] = VFP1_CM[2];
2597 IEV_TAB = IEV_TAB + 1;
2598
2599 if (AFP2 > 0)
2600 {
2601 EV_TAB[IEV_TAB][0] = ZFP2;
2602 EV_TAB[IEV_TAB][1] = AFP2;
2603 EV_TAB[IEV_TAB][5] = SFP2;
2604 EV_TAB[IEV_TAB][2] = VFP2_CM[0];
2605 EV_TAB[IEV_TAB][3] = VFP2_CM[1];
2606 EV_TAB[IEV_TAB][4] = VFP2_CM[2];
2607 IEV_TAB = IEV_TAB + 1;
2608 }
2609
2610 if (AFPIMF > 0)
2611 {
2612 EV_TAB[IEV_TAB][0] = ZFPIMF;
2613 EV_TAB[IEV_TAB][1] = AFPIMF;
2614 EV_TAB[IEV_TAB][5] = SFPIMF;
2615 EV_TAB[IEV_TAB][2] = VIMF_CM[0];
2616 EV_TAB[IEV_TAB][3] = VIMF_CM[1];
2617 EV_TAB[IEV_TAB][4] = VIMF_CM[2];
2618 IEV_TAB = IEV_TAB + 1;
2619 }
2620
2622 return;
2623}
constexpr const G4int indexpart
G4double C(G4double temp)
G4double eflmac(G4int ia, G4int iz, G4int flag, G4int optshp)
void evapora(G4double zprf, G4double aprf, G4double *ee_par, G4double jprf, G4double *zf_par, G4double *af_par, G4double *mtota_par, G4double *vleva_par, G4double *vxeva_par, G4double *vyeva_par, G4int *ff_par, G4int *fimf_par, G4double *fzimf, G4double *faimf, G4double *tkeimf_par, G4double *jprfout, G4int *inttype_par, G4int *inum_par, G4double EV_TEMP[indexpart][6], G4int *iev_tab_temp_par, G4int *nblam0)
G4int IPOWERLIMHAZ(G4double lambda, G4int xmin, G4int xmax)
void isostab_lim(G4int z, G4int *nmin, G4int *nmax)
G4double DSIGN(G4double a, G4double b)
G4int ISIGN(G4int a, G4int b)
void AMOMENT(G4double AABRA, G4double APRF, G4int IMULTIFR, G4double *PX, G4double *PY, G4double *PZ)
void tke_bu(G4double Z, G4double A, G4double ZALL, G4double AAL, G4double *VX, G4double *VY, G4double *VZ)
G4double dint(G4double a)
void fission(G4double AF, G4double ZF, G4double EE, G4double JPRF, G4double *VX1_FISSION, G4double *VY1_FISSION, G4double *VZ1_FISSION, G4double *VX2_FISSION, G4double *VY2_FISSION, G4double *VZ2_FISSION, G4int *ZFP1, G4int *AFP1, G4int *SFP1, G4int *ZFP2, G4int *AFP2, G4int *SFP2, G4int *imode, G4double *VX_EVA_SC, G4double *VY_EVA_SC, G4double *VZ_EVA_SC, G4double EV_TEMP[indexpart][6], G4int *IEV_TAB_FIS, G4int *NbLam0)
G4int max(G4int a, G4int b)
void unstable_nuclei(G4int AFP, G4int ZFP, G4int *AFPNEW, G4int *ZFPNEW, G4int &IOUNSTABLE, G4double VX, G4double VY, G4double VZ, G4double *VP1X, G4double *VP1Y, G4double *VP1Z, G4double BU_TAB_TEMP[indexpart][6], G4int *ILOOP)
void FillData(G4int IMULTBU, G4int IEV_TAB)
void SetParametersG4(G4int z, G4int a)
void lorentz_boost(G4double VXRIN, G4double VYRIN, G4double VZRIN, G4double VXIN, G4double VYIN, G4double VZIN, G4double *VXOUT, G4double *VYOUT, G4double *VZOUT)
1.13.2