Main interface to the de-excitation code for hyper-nuclei.
88{
91
93
94mult10:
96
97 varntp->clear();
98
99 if (nucleusS > 0) nucleusS = 0;
100
101 G4int NbLam0 = std::abs(nucleusS);
102
103 Ainit = -1 * nucleusA;
104 Zinit = -1 * nucleusZ;
105 Sinit = -1 * nucleusS;
106
109 G4int ZFP1 = 0, AFP1 = 0, AFPIMF = 0, ZFPIMF = 0, ZFP2 = 0, AFP2 = 0, SFP1 = 0, SFP2 = 0,
110 SFPIMF = 0;
111 G4double vx_eva = 0.0, vy_eva = 0.0, vz_eva = 0.0;
112 G4double VX_PREF = 0., VY_PREF = 0., VZ_PREF = 00, VP1X, VP1Y, VP1Z, VXOUT, VYOUT, VZOUT, V_CM[3],
113 VFP1_CM[3], VFP2_CM[3], VIMF_CM[3], VX2OUT, VY2OUT, VZ2OUT;
114 G4double zf = 0.0, af = 0.0, mtota = 0.0, tkeimf = 0.0, jprf0 = 0.;
115 G4int ff = 0, afpnew = 0, zfpnew = 0, aprfp = 0, zprfp = 0, IOUNSTABLE = 0, ILOOP = 0,
116 IEV_TAB = 0, IEV_TAB_TEMP = 0;
117 G4int fimf = 0, INMIN = 0, INMAX = 0;
119 G4int inum = eventnumber;
121 opt->optimfallowed = 1;
122
123 if (fiss->zt > 56) {
124 fiss->ifis = 1;
125 }
126 else {
127 fiss->ifis = 0;
128 }
129
130 if (NbLam0 > 0) {
131 opt->nblan0 = NbLam0;
132 }
133
138
143
144 gammaemission = 0;
145 G4double T_init = 0., T_diff = 0., a_tilda = 0., a_tilda_BU = 0., EE_diff = 0., EINCL = 0.,
146 A_FINAL = 0., Z_FINAL = 0., E_FINAL = 0.;
147
148 G4double A_diff = 0., ASLOPE1, ASLOPE2, A_ACC, ABU_SLOPE, ABU_SUM = 0., AMEM = 0., ZMEM = 0.,
149 EMEM = 0., JMEM = 0., PX_BU_SUM = 0.0, PY_BU_SUM = 0.0, PZ_BU_SUM = 0.0, ETOT_SUM = 0.,
150 P_BU_SUM = 0., ZBU_SUM = 0., Z_Breakup_sum = 0., A_Breakup, Z_Breakup, N_Breakup, G_SYMM,
151 CZ, Sigma_Z, Z_Breakup_Mean, ZTEMP = 0., ATEMP = 0.;
152
153 G4double ETOT_PRF = 0.0, PXPRFP = 0., PYPRFP = 0., PZPRFP = 0., PPRFP = 0., VX1_BU = 0.,
154 VY1_BU = 0., VZ1_BU = 0., VBU2 = 0., GAMMA_REL = 1.0, Eexc_BU_SUM = 0., VX_BU_SUM = 0.,
155 VY_BU_SUM = 0., VZ_BU_SUM = 0., E_tot_BU = 0., EKIN_BU = 0., ZIMFBU = 0., AIMFBU = 0.,
156 ZFFBU = 0., AFFBU = 0., AFBU = 0., ZFBU = 0., EEBU = 0., TKEIMFBU = 0., vx_evabu = 0.,
157 vy_evabu = 0., vz_evabu = 0., Bvalue_BU = 0., P_BU = 0., ETOT_BU = 1., PX_BU = 0.,
158 PY_BU = 0., PZ_BU = 0., VX2_BU = 0., VY2_BU = 0., VZ2_BU = 0.;
159
160 G4int ABU_DIFF, ZBU_DIFF, NBU_DIFF;
161 G4int INEWLOOP = 0, ILOOPBU = 0;
162
166
167 if (nucleusA < 1) {
168 std::cout << "Error - Remnant with a mass number A below 1." << std::endl;
169
170 return;
171 }
172
173 for (
G4int j = 0; j < 3; j++) {
174 V_CM[j] = 0.;
175 VFP1_CM[j] = 0.;
176 VFP2_CM[j] = 0.;
177 VIMF_CM[j] = 0.;
178 }
179
181 for (
G4int I2 = 0; I2 < 12; I2++)
182 BU_TAB[I1][I2] = 0.0;
183 for (
G4int I2 = 0; I2 < 6; I2++) {
184 BU_TAB_TEMP[I1][I2] = 0.0;
185 BU_TAB_TEMP1[I1][I2] = 0.0;
186 EV_TAB_TEMP[I1][I2] = 0.0;
187 EV_TAB[I1][I2] = 0.0;
188 EV_TAB_SSC[I1][I2] = 0.0;
189 EV_TEMP[I1][I2] = 0.0;
190 }
191 }
192
194 if (idebug == 1) {
195 zprf = 81.;
196 aprf = 201.;
197
198 ee = 100.0;
199 jprf = 10.;
200 zf = 0.;
201 af = 0.;
202 mtota = 0.;
203 ff = 1;
204 inttype = 0;
205
206 }
207
210 EINCL = ee;
211
212
213
214
215
216 G4double pincl = std::sqrt(pxrem * pxrem + pyrem * pyrem + pzrem * pzrem);
217
218 G4double ETOT_incl = std::sqrt(pincl * pincl + (AAINCL * amu) * (AAINCL * amu));
219 G4double VX_incl =
C * pxrem / ETOT_incl;
220 G4double VY_incl =
C * pyrem / ETOT_incl;
221 G4double VZ_incl =
C * pzrem / ETOT_incl;
222
223
228 IEV_TAB = 0;
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252 T_freeze_out =
253 (T_freeze_out_in >= 0.0) ? T_freeze_out_in :
max(9.33 * std::exp(-0.00282 * AAINCL), 5.5);
254
255 a_tilda =
256 ald->av * aprf + ald->as * std::pow(aprf, 2.0 / 3.0) + ald->ak * std::pow(aprf, 1.0 / 3.0);
257
258 T_init = std::sqrt(EINCL / a_tilda);
259
260 T_diff = T_init - T_freeze_out;
261
262 if (T_diff > 0.1 && zprf > 2. && (aprf - zprf) > 0.) {
263
264
265 varntp->kfis = 10;
266
267 for (
G4int i = 0; i < 5; i++) {
268 EE_diff = EINCL - a_tilda * T_freeze_out * T_freeze_out;
269
270
271
272
273
274
275 A_diff =
dint(EE_diff / (8.0 * 5.0 / T_freeze_out));
276
277 if (A_diff > AAINCL) A_diff = AAINCL;
278
279 A_FINAL = AAINCL - A_diff;
280
281 a_tilda = ald->av * A_FINAL + ald->as * std::pow(A_FINAL, 2.0 / 3.0)
282 + ald->ak * std::pow(A_FINAL, 1.0 / 3.0);
283 E_FINAL = a_tilda * T_freeze_out * T_freeze_out;
284
285 if (A_FINAL < 4.0) {
286 EE_diff = EINCL - E_FINAL;
287 A_FINAL = 1.0;
288 Z_FINAL = 1.0;
289 E_FINAL = 0.0;
290 goto mul4325;
291 }
292 }
293 mul4325:
294
295
296
297
298 Z_FINAL =
dint(zprf * A_FINAL / (aprf));
299
300 if (E_FINAL < 0.0) E_FINAL = 0.0;
301
302 aprf = A_FINAL;
303 zprf = Z_FINAL;
304 ee = E_FINAL;
305
306 A_diff = AAINCL - aprf;
307
308
309 if (A_diff <= 1.0) {
310 aprf = AAINCL;
311 zprf = ZAINCL;
312 ee = EINCL;
313 IMULTIFR = 0;
314 goto mult7777;
315 }
316 else if (A_diff > 1.0) {
317 A_ACC = 0.0;
318
319
320 ASLOPE1 = -2.400;
321 ASLOPE2 = -1.200;
322
323 a_tilda = ald->av * AAINCL + ald->as * std::pow(AAINCL, 2.0 / 3.0)
324 + ald->ak * std::pow(AAINCL, 1.0 / 3.0);
325
326 E_FINAL = a_tilda * T_freeze_out * T_freeze_out;
327
328 ABU_SLOPE =
329 (ASLOPE1 - ASLOPE2) / 4.0 * (E_FINAL / AAINCL) + ASLOPE1 - (ASLOPE1 - ASLOPE2) * 7.0 / 4.0;
330
331
332
333
334
335
336
337
338
339
340
341 if (ABU_SLOPE > -1.01) ABU_SLOPE = -1.01;
342
343 I_Breakup = 0;
344 Z_Breakup_sum = Z_FINAL;
345 ABU_SUM = 0.0;
346 ZBU_SUM = 0.0;
347
348 for (
G4int i = 0; i < 100; i++) {
349 IS = 0;
350 mult4326:
352
353 IS = IS + 1;
354 if (IS > 100) {
355 std::cout << "WARNING: IPOWERLIMHAZ CALLED MORE THAN 100 TIMES WHEN "
356 "CALCULATING A_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED: "
357 << A_Breakup << std::endl;
358 goto mult10;
359 }
360
361 if (A_Breakup > AAINCL) goto mult4326;
362
363 if (A_Breakup <= 0.0) {
364 std::cout << "A_BREAKUP <= 0 " << std::endl;
365 goto mult10;
366 }
367
368 A_ACC = A_ACC + A_Breakup;
369
370 if (A_ACC <= A_diff) {
371 Z_Breakup_Mean =
dint(A_Breakup * ZAINCL / AAINCL);
372
373 Z_Breakup_sum = Z_Breakup_sum + Z_Breakup_Mean;
374
375
376 G_SYMM = 34.2281 - 5.14037 * E_FINAL / AAINCL;
377 if (E_FINAL / AAINCL < 2.0) G_SYMM = 25.0;
378 if (E_FINAL / AAINCL > 4.0) G_SYMM = 15.0;
379
380
381
382 G_SYMM = 25.0;
383 CZ = 2.0 * G_SYMM * 4.0 / A_Breakup;
384
385
387 Sigma_Z = std::sqrt(T_freeze_out / CZ);
388
389 IS = 0;
390 mult4333:
392 IS = IS + 1;
393
394 if (IS > 100) {
395 std::cout << "WARNING: GAUSSHAZ CALLED MORE THAN 100 TIMES WHEN "
396 "CALCULATING Z_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE "
397 "DICED: "
398 << A_Breakup << " " << Z_Breakup << std::endl;
399 goto mult10;
400 }
401
402 if (Z_Breakup < 0.0) goto mult4333;
403 if ((A_Breakup - Z_Breakup) < 0.0) goto mult4333;
404 if ((A_Breakup - Z_Breakup) == 0.0 && Z_Breakup != 1.0) goto mult4333;
405
406 if (Z_Breakup >= ZAINCL) {
407 IIS = IIS + 1;
408 if (IIS > 10) {
409 std::cout << "Z_BREAKUP RESAMPLED MORE THAN 10 TIMES; EVENT WILL "
410 "BE RESAMPLED AGAIN "
411 << std::endl;
412 goto mult10;
413 }
414 goto mult4333;
415 }
416
417
419
420 if (Z_Breakup > 2.0) {
421 if (
idnint(A_Breakup - Z_Breakup) < INMIN
422 ||
idnint(A_Breakup - Z_Breakup) > (INMAX + 5))
423 {
424
425
426 goto mult4343;
427 }
428 }
429
430 mult4343:
431
432
433
434
435 N_Breakup = A_Breakup - Z_Breakup;
436 BU_TAB[I_Breakup][0] =
dint(Z_Breakup);
437 BU_TAB[I_Breakup][1] =
dint(A_Breakup);
438 ABU_SUM = ABU_SUM + BU_TAB[i][1];
439 ZBU_SUM = ZBU_SUM + BU_TAB[i][0];
440
441
442 BU_TAB[I_Breakup][3] = 0.0;
443 I_Breakup = I_Breakup + 1;
444 IMULTBU = IMULTBU + 1;
445 }
446 else {
447
448
449
450
451
452
453 goto mult4327;
454 }
455 }
456
457
458 }
459 mult4327:
460 IMULTIFR = 1;
461
462
463 ABU_DIFF =
idnint(ABU_SUM + aprf - AAINCL);
464 ZBU_DIFF =
idnint(ZBU_SUM + zprf - ZAINCL);
465 NBU_DIFF =
idnint((ABU_SUM - ZBU_SUM) + (aprf - zprf) - (AAINCL - ZAINCL));
466
467 if (IMULTBU > 200) std::cout << "WARNING - MORE THAN 200 BU " << IMULTBU << std::endl;
468
469 if (IMULTBU < 1) std::cout << "WARNING - LESS THAN 1 BU " << IMULTBU << std::endl;
470
471
473 for (
G4int i = 0; i < IMULTBU; i++)
474 IMEM_BU[i] = 0;
475
476 while (NBU_DIFF != 0 && ZBU_DIFF != 0) {
477
478
479
480 IS = 0;
481 mult5555:
483 IPROBA = IPROBA + 1;
484 IS = IS + 1;
485 if (IS > 100) {
486 std::cout << "WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING "
487 "N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED."
488 << std::endl;
489 goto mult10;
490 }
492 if (IMEM_BU[IEL] == 1) goto mult5555;
493 if (!(IEL < 200)) std::cout << "5555:" << IEL << RHAZ << IMULTBU << std::endl;
494 if (IEL < 0) std::cout << "5555:" << IEL << RHAZ << IMULTBU << std::endl;
495 if (IEL <= IMULTBU) {
496 N_Breakup =
dint(BU_TAB[IEL][1] - BU_TAB[IEL][0] -
DSIGN(1.0,
G4double(NBU_DIFF)));
497 }
498 else if (IEL > IMULTBU) {
500 }
501 if (N_Breakup < 0.0) {
502 IMEM_BU[IEL] = 1;
503 goto mult5555;
504 }
505 if (IEL <= IMULTBU) {
507 }
508 else if (IEL > IMULTBU) {
510 }
511 if (ZTEMP < 0.0) {
512 IMEM_BU[IEL] = 1;
513 goto mult5555;
514 }
515 if (ZTEMP < 1.0 && N_Breakup < 1.0) {
516 IMEM_BU[IEL] = 1;
517 goto mult5555;
518 }
519
520
521
522
523
524
525
526 if (IEL <= IMULTBU) {
527 BU_TAB[IEL][0] =
dint(ZTEMP);
528 BU_TAB[IEL][1] =
dint(ZTEMP + N_Breakup);
529 }
530 else if (IEL > IMULTBU) {
532 aprf =
dint(ZTEMP + N_Breakup);
533 }
534 NBU_DIFF = NBU_DIFF -
ISIGN(1, NBU_DIFF);
535 ZBU_DIFF = ZBU_DIFF -
ISIGN(1, ZBU_DIFF);
536 }
537
538 IPROBA = 0;
539 for (
G4int i = 0; i < IMULTBU; i++)
540 IMEM_BU[i] = 0;
541
542 if (NBU_DIFF != 0 && ZBU_DIFF == 0) {
543 while (NBU_DIFF > 0 || NBU_DIFF < 0) {
544 IS = 0;
545 mult5556:
547 IS = IS + 1;
548 if (IS > 100) {
549 std::cout << "WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING "
550 "N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED."
551 << std::endl;
552 goto mult10;
553 }
555 if (IMEM_BU[IEL] == 1) goto mult5556;
556
557 if (IPROBA > IMULTBU + 1 && NBU_DIFF > 0) {
558 std::cout << "###',IPROBA,IMULTBU,NBU_DIFF,ZBU_DIFF,T_freeze_out" << std::endl;
559 IPROBA = IPROBA + 1;
560 if (IEL <= IMULTBU) {
561 BU_TAB[IEL][1] =
dint(BU_TAB[IEL][1] -
G4double(NBU_DIFF));
562 }
563 else {
564 if (IEL > IMULTBU) aprf =
dint(aprf -
G4double(NBU_DIFF));
565 }
566 goto mult5432;
567 }
568 if (!(IEL < 200)) std::cout << "5556:" << IEL << RHAZ << IMULTBU << std::endl;
569 if (IEL < 0) std::cout << "5556:" << IEL << RHAZ << IMULTBU << std::endl;
570 if (IEL <= IMULTBU) {
571 N_Breakup =
dint(BU_TAB[IEL][1] - BU_TAB[IEL][0] -
DSIGN(1.0,
G4double(NBU_DIFF)));
572 }
573 else if (IEL > IMULTBU) {
575 }
576 if (N_Breakup < 0.0) {
577 IMEM_BU[IEL] = 1;
578 goto mult5556;
579 }
580 if (IEL <= IMULTBU) {
581 ATEMP =
dint(BU_TAB[IEL][0] + N_Breakup);
582 }
583 else if (IEL > IMULTBU) {
584 ATEMP =
dint(zprf + N_Breakup);
585 }
586 if ((ATEMP - N_Breakup) < 1.0 && N_Breakup < 1.0) {
587 IMEM_BU[IEL] = 1;
588 goto mult5556;
589 }
590
591
592
593
594
595 if (IEL <= IMULTBU)
596 BU_TAB[IEL][1] =
dint(BU_TAB[IEL][0] + N_Breakup);
597 else if (IEL > IMULTBU)
598 aprf =
dint(zprf + N_Breakup);
599
600 NBU_DIFF = NBU_DIFF -
ISIGN(1, NBU_DIFF);
601 }
602
603 IPROBA = 0;
604 for (
G4int i = 0; i < IMULTBU; i++)
605 IMEM_BU[i] = 0;
606 }
607 else {
608 if (ZBU_DIFF != 0 && NBU_DIFF == 0) {
609 while (ZBU_DIFF > 0 || ZBU_DIFF < 0) {
610 IS = 0;
611 mult5557:
613 IS = IS + 1;
614 if (IS > 100) {
615 std::cout << "WARNING: HAZ CALLED MORE THAN 100 TIMES WHEN CALCULATING "
616 "N_BREAKUP IN Rn07.FOR. NEW EVENT WILL BE DICED."
617 << std::endl;
618 goto mult10;
619 }
621 if (IMEM_BU[IEL] == 1) goto mult5557;
622
623 if (IPROBA > IMULTBU + 1 && ZBU_DIFF > 0) {
624 std::cout << "###',IPROBA,IMULTBU,NBU_DIFF,ZBU_DIFF,T_freeze_out" << std::endl;
625 IPROBA = IPROBA + 1;
626 if (IEL <= IMULTBU) {
627 N_Breakup =
dint(BU_TAB[IEL][1] - BU_TAB[IEL][0]);
628 BU_TAB[IEL][0] =
dint(BU_TAB[IEL][0] -
G4double(ZBU_DIFF));
629 BU_TAB[IEL][1] =
dint(BU_TAB[IEL][0] + N_Breakup);
630 }
631 else {
632 if (IEL > IMULTBU) {
633 N_Breakup = aprf - zprf;
635 aprf =
dint(zprf + N_Breakup);
636 }
637 }
638 goto mult5432;
639 }
640 if (!(IEL < 200)) std::cout << "5557:" << IEL << RHAZ << IMULTBU << std::endl;
641 if (IEL < 0) std::cout << "5557:" << IEL << RHAZ << IMULTBU << std::endl;
642 if (IEL <= IMULTBU) {
643 N_Breakup =
dint(BU_TAB[IEL][1] - BU_TAB[IEL][0]);
645 }
646 else if (IEL > IMULTBU) {
647 N_Breakup =
dint(aprf - zprf);
649 }
650 ATEMP =
dint(ZTEMP + N_Breakup);
651 if (ZTEMP < 0.0) {
652 IMEM_BU[IEL] = 1;
653 goto mult5557;
654 }
655 if ((ATEMP - ZTEMP) < 0.0) {
656 IMEM_BU[IEL] = 1;
657 goto mult5557;
658 }
659 if ((ATEMP - ZTEMP) < 1.0 && ZTEMP < 1.0) {
660 IMEM_BU[IEL] = 1;
661 goto mult5557;
662 }
663 if (IEL <= IMULTBU) {
664 BU_TAB[IEL][0] =
dint(ZTEMP);
665 BU_TAB[IEL][1] =
dint(ZTEMP + N_Breakup);
666 }
667 else {
668 if (IEL > IMULTBU) {
670 aprf =
dint(ZTEMP + N_Breakup);
671 }
672 }
673 ZBU_DIFF = ZBU_DIFF -
ISIGN(1, ZBU_DIFF);
674 }
675 }
676 }
677
678 mult5432:
679
680
681 ZMEM = 0.0;
682
683 for (
G4int i = 0; i < IMULTBU; i++) {
684
685
686
687
688
689 if (BU_TAB[i][0] > 2.0) {
690 a_tilda_BU = ald->av * BU_TAB[i][1] + ald->as * std::pow(BU_TAB[i][1], 2.0 / 3.0)
691 + ald->ak * std::pow(BU_TAB[i][1], 1.0 / 3.0);
692 BU_TAB[i][2] = a_tilda_BU * T_freeze_out * T_freeze_out;
693 }
694 else {
695 BU_TAB[i][2] = 0.0;
696 }
697
698 if (BU_TAB[i][0] > ZMEM) {
699 IMEM = i;
700 ZMEM = BU_TAB[i][0];
701 AMEM = BU_TAB[i][1];
702 EMEM = BU_TAB[i][2];
703 JMEM = BU_TAB[i][3];
704 }
705 }
706
707 if (zprf < ZMEM) {
708 BU_TAB[IMEM][0] = zprf;
709 BU_TAB[IMEM][1] = aprf;
710 BU_TAB[IMEM][2] = ee;
711 BU_TAB[IMEM][3] = jprf;
712 zprf = ZMEM;
713 aprf = AMEM;
716 ee = EMEM;
717 jprf = JMEM;
718 }
719
720
721 ABU_SUM = aprf;
722 ZBU_SUM = zprf;
723 for (
G4int i = 0; i < IMULTBU; i++) {
724 ABU_SUM = ABU_SUM + BU_TAB[i][1];
725 ZBU_SUM = ZBU_SUM + BU_TAB[i][0];
726 }
727 ABU_DIFF =
idnint(ABU_SUM - AAINCL);
728 ZBU_DIFF =
idnint(ZBU_SUM - ZAINCL);
729
730 if (ABU_DIFF != 0 || ZBU_DIFF != 0)
731 std::cout << "Problem of mass in BU " << ABU_DIFF << " " << ZBU_DIFF << std::endl;
732 PX_BU_SUM = 0.0;
733 PY_BU_SUM = 0.0;
734 PZ_BU_SUM = 0.0;
735
736
737
738 AMOMENT(AAINCL, aprf, 1, &PXPRFP, &PYPRFP, &PZPRFP);
739 PPRFP = std::sqrt(PXPRFP * PXPRFP + PYPRFP * PYPRFP + PZPRFP * PZPRFP);
740
741
742 ETOT_PRF = std::sqrt(PPRFP * PPRFP + (aprf * amu) * (aprf * amu));
743 VX_PREF =
C * PXPRFP / ETOT_PRF;
744 VY_PREF =
C * PYPRFP / ETOT_PRF;
745 VZ_PREF =
C * PZPRFP / ETOT_PRF;
746
747
748 tke_bu(zprf, aprf, ZAINCL, AAINCL, &VX1_BU, &VY1_BU, &VZ1_BU);
749
750
751
752
753
754
755 lorentz_boost(VX1_BU, VY1_BU, VZ1_BU, VX_PREF, VY_PREF, VZ_PREF, &VXOUT, &VYOUT, &VZOUT);
756
757 VX_PREF = VXOUT;
758 VY_PREF = VYOUT;
759 VZ_PREF = VZOUT;
760
761
762 VBU2 = VX_PREF * VX_PREF + VY_PREF * VY_PREF + VZ_PREF * VZ_PREF;
763 GAMMA_REL = std::sqrt(1.0 - VBU2 / (
C *
C));
764 ETOT_PRF = aprf * amu / GAMMA_REL;
765 PXPRFP = ETOT_PRF * VX_PREF /
C;
766 PYPRFP = ETOT_PRF * VY_PREF /
C;
767 PZPRFP = ETOT_PRF * VZ_PREF /
C;
768
769
770
771
772
773
774 PX_BU_SUM = PXPRFP;
775 PY_BU_SUM = PYPRFP;
776 PZ_BU_SUM = PZPRFP;
777
778 Eexc_BU_SUM = ee;
780
781 for (I_Breakup = 0; I_Breakup < IMULTBU; I_Breakup++) {
782
783 Bvalue_BU =
785 Eexc_BU_SUM = Eexc_BU_SUM + BU_TAB[I_Breakup][2];
786
787 AMOMENT(AAINCL, BU_TAB[I_Breakup][1], 1, &PX_BU, &PY_BU, &PZ_BU);
788 P_BU = std::sqrt(PX_BU * PX_BU + PY_BU * PY_BU + PZ_BU * PZ_BU);
789
790
791 ETOT_BU =
792 std::sqrt(P_BU * P_BU + (BU_TAB[I_Breakup][1] * amu) * (BU_TAB[I_Breakup][1] * amu));
793 BU_TAB[I_Breakup][4] =
C * PX_BU / ETOT_BU;
794 BU_TAB[I_Breakup][5] =
C * PY_BU / ETOT_BU;
795 BU_TAB[I_Breakup][6] =
C * PZ_BU / ETOT_BU;
796
797 tke_bu(BU_TAB[I_Breakup][0], BU_TAB[I_Breakup][1], ZAINCL, AAINCL, &VX2_BU, &VY2_BU, &VZ2_BU);
798
799
800
801
802
803
804 lorentz_boost(VX2_BU, VY2_BU, VZ2_BU, BU_TAB[I_Breakup][4], BU_TAB[I_Breakup][5],
805 BU_TAB[I_Breakup][6], &VXOUT, &VYOUT, &VZOUT);
806
807 BU_TAB[I_Breakup][4] = VXOUT;
808 BU_TAB[I_Breakup][5] = VYOUT;
809 BU_TAB[I_Breakup][6] = VZOUT;
810
811
812 VBU2 = BU_TAB[I_Breakup][4] * BU_TAB[I_Breakup][4]
813 + BU_TAB[I_Breakup][5] * BU_TAB[I_Breakup][5]
814 + BU_TAB[I_Breakup][6] * BU_TAB[I_Breakup][6];
815 GAMMA_REL = std::sqrt(1.0 - VBU2 / (
C *
C));
816 ETOT_BU = BU_TAB[I_Breakup][1] * amu / GAMMA_REL;
817 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
818 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
819 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
820
821 PX_BU_SUM = PX_BU_SUM + PX_BU;
822 PY_BU_SUM = PY_BU_SUM + PY_BU;
823 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
824
825 }
826
827
828 P_BU_SUM = std::sqrt(PX_BU_SUM * PX_BU_SUM + PY_BU_SUM * PY_BU_SUM + PZ_BU_SUM * PZ_BU_SUM);
829
830
831 ETOT_SUM = std::sqrt(P_BU_SUM * P_BU_SUM + (AAINCL * amu) * (AAINCL * amu));
832
833 VX_BU_SUM =
C * PX_BU_SUM / ETOT_SUM;
834 VY_BU_SUM =
C * PY_BU_SUM / ETOT_SUM;
835 VZ_BU_SUM =
C * PZ_BU_SUM / ETOT_SUM;
836
837
838
839
840
841
842 lorentz_boost(-VX_BU_SUM, -VY_BU_SUM, -VZ_BU_SUM, VX_PREF, VY_PREF, VZ_PREF, &VXOUT, &VYOUT,
843 &VZOUT);
844
845 VX_PREF = VXOUT;
846 VY_PREF = VYOUT;
847 VZ_PREF = VZOUT;
848
849 VBU2 = VX_PREF * VX_PREF + VY_PREF * VY_PREF + VZ_PREF * VZ_PREF;
850 GAMMA_REL = std::sqrt(1.0 - VBU2 / (
C *
C));
851 ETOT_PRF = aprf * amu / GAMMA_REL;
852 PXPRFP = ETOT_PRF * VX_PREF /
C;
853 PYPRFP = ETOT_PRF * VY_PREF /
C;
854 PZPRFP = ETOT_PRF * VZ_PREF /
C;
855
856 PX_BU_SUM = 0.0;
857 PY_BU_SUM = 0.0;
858 PZ_BU_SUM = 0.0;
859
860 PX_BU_SUM = PXPRFP;
861 PY_BU_SUM = PYPRFP;
862 PZ_BU_SUM = PZPRFP;
863 E_tot_BU = ETOT_PRF;
864
865 EKIN_BU = aprf * amu / GAMMA_REL - aprf * amu;
866
867 for (I_Breakup = 0; I_Breakup < IMULTBU; I_Breakup++) {
868
869
870
871
872
873 lorentz_boost(-VX_BU_SUM, -VY_BU_SUM, -VZ_BU_SUM, BU_TAB[I_Breakup][4], BU_TAB[I_Breakup][5],
874 BU_TAB[I_Breakup][6], &VXOUT, &VYOUT, &VZOUT);
875
876 BU_TAB[I_Breakup][4] = VXOUT;
877 BU_TAB[I_Breakup][5] = VYOUT;
878 BU_TAB[I_Breakup][6] = VZOUT;
879
880 VBU2 = BU_TAB[I_Breakup][4] * BU_TAB[I_Breakup][4]
881 + BU_TAB[I_Breakup][5] * BU_TAB[I_Breakup][5]
882 + BU_TAB[I_Breakup][6] * BU_TAB[I_Breakup][6];
883 GAMMA_REL = std::sqrt(1.0 - VBU2 / (
C *
C));
884
885 ETOT_BU = BU_TAB[I_Breakup][1] * amu / GAMMA_REL;
886
887 EKIN_BU = EKIN_BU + BU_TAB[I_Breakup][1] * amu / GAMMA_REL - BU_TAB[I_Breakup][1] * amu;
888
889 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
890 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
891 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
892 E_tot_BU = E_tot_BU + ETOT_BU;
893
894 PX_BU_SUM = PX_BU_SUM + PX_BU;
895 PY_BU_SUM = PY_BU_SUM + PY_BU;
896 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
897 }
898
899 if (std::abs(PX_BU_SUM) > 10. || std::abs(PY_BU_SUM) > 10. || std::abs(PZ_BU_SUM) > 10.) {
900
901 P_BU_SUM = std::sqrt(PX_BU_SUM * PX_BU_SUM + PY_BU_SUM * PY_BU_SUM + PZ_BU_SUM * PZ_BU_SUM);
902
903
904 ETOT_SUM = std::sqrt(P_BU_SUM * P_BU_SUM + (AAINCL * amu) * (AAINCL * amu));
905
906 VX_BU_SUM =
C * PX_BU_SUM / ETOT_SUM;
907 VY_BU_SUM =
C * PY_BU_SUM / ETOT_SUM;
908 VZ_BU_SUM =
C * PZ_BU_SUM / ETOT_SUM;
909
910
911
912
913
914
915 lorentz_boost(-VX_BU_SUM, -VY_BU_SUM, -VZ_BU_SUM, VX_PREF, VY_PREF, VZ_PREF, &VXOUT, &VYOUT,
916 &VZOUT);
917
918 VX_PREF = VXOUT;
919 VY_PREF = VYOUT;
920 VZ_PREF = VZOUT;
921
922 VBU2 = VX_PREF * VX_PREF + VY_PREF * VY_PREF + VZ_PREF * VZ_PREF;
923 GAMMA_REL = std::sqrt(1.0 - VBU2 / (
C *
C));
924 ETOT_PRF = aprf * amu / GAMMA_REL;
925 PXPRFP = ETOT_PRF * VX_PREF /
C;
926 PYPRFP = ETOT_PRF * VY_PREF /
C;
927 PZPRFP = ETOT_PRF * VZ_PREF /
C;
928
929 PX_BU_SUM = 0.0;
930 PY_BU_SUM = 0.0;
931 PZ_BU_SUM = 0.0;
932
933 PX_BU_SUM = PXPRFP;
934 PY_BU_SUM = PYPRFP;
935 PZ_BU_SUM = PZPRFP;
936 E_tot_BU = ETOT_PRF;
937
938 EKIN_BU = aprf * amu / GAMMA_REL - aprf * amu;
939
940 for (I_Breakup = 0; I_Breakup < IMULTBU; I_Breakup++) {
941
942
943
944
945
946 lorentz_boost(-VX_BU_SUM, -VY_BU_SUM, -VZ_BU_SUM, BU_TAB[I_Breakup][4],
947 BU_TAB[I_Breakup][5], BU_TAB[I_Breakup][6], &VXOUT, &VYOUT, &VZOUT);
948
949 BU_TAB[I_Breakup][4] = VXOUT;
950 BU_TAB[I_Breakup][5] = VYOUT;
951 BU_TAB[I_Breakup][6] = VZOUT;
952
953 VBU2 = BU_TAB[I_Breakup][4] * BU_TAB[I_Breakup][4]
954 + BU_TAB[I_Breakup][5] * BU_TAB[I_Breakup][5]
955 + BU_TAB[I_Breakup][6] * BU_TAB[I_Breakup][6];
956 GAMMA_REL = std::sqrt(1.0 - VBU2 / (
C *
C));
957
958 ETOT_BU = BU_TAB[I_Breakup][1] * amu / GAMMA_REL;
959
960 EKIN_BU = EKIN_BU + BU_TAB[I_Breakup][1] * amu / GAMMA_REL - BU_TAB[I_Breakup][1] * amu;
961
962 PX_BU = ETOT_BU * BU_TAB[I_Breakup][4] /
C;
963 PY_BU = ETOT_BU * BU_TAB[I_Breakup][5] /
C;
964 PZ_BU = ETOT_BU * BU_TAB[I_Breakup][6] /
C;
965 E_tot_BU = E_tot_BU + ETOT_BU;
966
967 PX_BU_SUM = PX_BU_SUM + PX_BU;
968 PY_BU_SUM = PY_BU_SUM + PY_BU;
969 PZ_BU_SUM = PZ_BU_SUM + PZ_BU;
970 }
971 }
972
973
974
975
976
977 INEWLOOP = 0;
978 for (
G4int i = 0; i < IMULTBU; i++) {
979 if (BU_TAB[i][0] < 3.0 || BU_TAB[i][0] == BU_TAB[i][1]) {
981 BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VP1X, &VP1Y, &VP1Z, BU_TAB_TEMP,
982 &ILOOP);
983
984 if (IOUNSTABLE > 0) {
985
988 BU_TAB[i][4] = VP1X;
989 BU_TAB[i][5] = VP1Y;
990 BU_TAB[i][6] = VP1Z;
991
992
993 for (int IJ = 0; IJ < ILOOP; IJ++) {
994 BU_TAB[IMULTBU + INEWLOOP + IJ][0] = BU_TAB_TEMP[IJ][0];
995 BU_TAB[IMULTBU + INEWLOOP + IJ][1] = BU_TAB_TEMP[IJ][1];
996 BU_TAB[IMULTBU + INEWLOOP + IJ][4] = BU_TAB_TEMP[IJ][2];
997 BU_TAB[IMULTBU + INEWLOOP + IJ][5] = BU_TAB_TEMP[IJ][3];
998 BU_TAB[IMULTBU + INEWLOOP + IJ][6] = BU_TAB_TEMP[IJ][4];
999 BU_TAB[IMULTBU + INEWLOOP + IJ][2] = 0.0;
1000 BU_TAB[IMULTBU + INEWLOOP + IJ][3] = 0.0;
1001 }
1002
1003 INEWLOOP = INEWLOOP + ILOOP;
1004
1005 }
1006 }
1007 }
1008
1009
1010 IMULTBU = IMULTBU + INEWLOOP;
1011
1012 opt->optimfallowed = 1;
1013 fiss->ifis = 0;
1014 gammaemission = 0;
1015 ILOOPBU = 0;
1016
1017
1018 std::vector<G4double> problamb(IMULTBU);
1019
1021 for (
G4int i = 0; i < IMULTBU; i++)
1022 sumN = sumN + BU_TAB[i][1] - BU_TAB[i][0];
1023
1024 for (
G4int i = 0; i < IMULTBU; i++) {
1025 problamb[i] = (BU_TAB[i][1] - BU_TAB[i][0]) / sumN;
1026 }
1027
1028 std::vector<G4int> Nblamb(IMULTBU, 0);
1029 for (
G4int j = 0; j < NbLam0;) {
1030 G4double probtotal = (aprf - zprf) / sumN;
1032
1033 if (ran <= probtotal) {
1034 NbLamprf++;
1035 goto directlamb0;
1036 }
1037 for (
G4int i = 0; i < IMULTBU; i++) {
1038
1039 if (probtotal < ran && ran <= probtotal + problamb[i]) {
1040 Nblamb[i] = Nblamb[i] + 1;
1041 goto directlamb0;
1042 }
1043 probtotal = probtotal + problamb[i];
1044 }
1045 directlamb0:
1046 j++;
1047 }
1048
1049 for (
G4int i = 0; i < IMULTBU; i++) {
1050 EEBU = BU_TAB[i][2];
1051 BU_TAB[i][10] = BU_TAB[i][6];
1053 if (BU_TAB[i][0] > 2.0) {
1054 G4int nbl = Nblamb[i];
1055 evapora(BU_TAB[i][0], BU_TAB[i][1], &EEBU, 0.0, &ZFBU, &AFBU, &mtota, &vz_evabu, &vx_evabu,
1056 &vy_evabu, &ff, &fimf, &ZIMFBU, &AIMFBU, &TKEIMFBU, &jprfbu, &inttype, &inum,
1057 EV_TEMP, &IEV_TAB_TEMP, &nbl);
1058
1059 Nblamb[i] = nbl;
1060 BU_TAB[i][9] = jprfbu;
1061
1062
1063
1064 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++) {
1065 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1066 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1067 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1068
1069
1070
1071
1072
1073 lorentz_boost(BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], EV_TEMP[IJ][2], EV_TEMP[IJ][3],
1074 EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1075 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1076 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1077 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1078 }
1079 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1080
1081
1082
1083
1084
1085
1086
1087 lorentz_boost(vx_evabu, vy_evabu, vz_evabu, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6],
1088 &VXOUT, &VYOUT, &VZOUT);
1089 BU_TAB[i][4] = VXOUT;
1090 BU_TAB[i][5] = VYOUT;
1091 BU_TAB[i][6] = VZOUT;
1092
1093 if (fimf == 0) {
1094 BU_TAB[i][7] =
dint(ZFBU);
1095 BU_TAB[i][8] =
dint(AFBU);
1096 BU_TAB[i][11] = nbl;
1097 }
1098
1099 if (fimf == 1) {
1100
1101
1102
1103
1106 opt->optimfallowed = 0;
1107 fiss->ifis = 0;
1108
1109 G4double EkinR1 = TKEIMFBU * AIMFBU / (AFBU + AIMFBU);
1110 G4double EkinR2 = TKEIMFBU * AFBU / (AFBU + AIMFBU);
1111 G4double V1 = std::sqrt(EkinR1 / AFBU) * 1.3887;
1112 G4double V2 = std::sqrt(EkinR2 / AIMFBU) * 1.3887;
1114 G4double VPERP1 = std::sqrt(
V1 *
V1 - VZ1_IMF * VZ1_IMF);
1116 G4double VX1_IMF = VPERP1 * std::sin(ALPHA1);
1117 G4double VY1_IMF = VPERP1 * std::cos(ALPHA1);
1121
1122 G4double EEIMFP = EEBU * AFBU / (AFBU + AIMFBU);
1123 G4double EEIMF = EEBU * AIMFBU / (AFBU + AIMFBU);
1124
1125
1126 G4double IINERTTOT = 0.40 * 931.490 * 1.160 * 1.160
1127 * (std::pow(AIMFBU, 5.0 / 3.0) + std::pow(AFBU, 5.0 / 3.0))
1128 + 931.490 * 1.160 * 1.160 * AIMFBU * AFBU / (AIMFBU + AFBU)
1129 * (std::pow(AIMFBU, 1. / 3.) + std::pow(AFBU, 1. / 3.))
1130 * (std::pow(AIMFBU, 1. / 3.) + std::pow(AFBU, 1. / 3.));
1131
1133 BU_TAB[i][9] * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(AFBU, 5.0 / 3.0) / IINERTTOT;
1135 BU_TAB[i][9] * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(AIMFBU, 5.0 / 3.0) / IINERTTOT;
1136
1137
1138
1139
1140
1141
1142 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VXOUT,
1143 &VYOUT, &VZOUT);
1144 BU_TAB[i][4] = VXOUT;
1145 BU_TAB[i][5] = VYOUT;
1146 BU_TAB[i][6] = VZOUT;
1147
1148 G4double vx1ev_imf = 0., vy1ev_imf = 0., vz1ev_imf = 0., zdummy = 0., adummy = 0.,
1149 tkedummy = 0., jprf1 = 0.;
1150
1151
1154 G4double pbH = (AFBU - ZFBU) / (AFBU - ZFBU + AIMFBU - ZIMFBU);
1155 for (
G4int j = 0; j < nbl; j++) {
1157 NbLamH++;
1158 }
1159 else {
1160 NbLamimf++;
1161 }
1162 }
1163
1164 evapora(ZFBU, AFBU, &EEIMFP, JPRFHEAVY, &ZFFBU, &AFFBU, &mtota, &vz1ev_imf, &vx1ev_imf,
1165 &vy1ev_imf, &FFBU1, &FIMFBU1, &zdummy, &adummy, &tkedummy, &jprf1, &inttype,
1166 &inum, EV_TEMP, &IEV_TAB_TEMP, &NbLamH);
1167
1168 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++) {
1169 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1170 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1171 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1172
1173
1174
1175
1176
1177 lorentz_boost(BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], EV_TEMP[IJ][2], EV_TEMP[IJ][3],
1178 EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1179 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1180 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1181 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1182 }
1183 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1184
1185 BU_TAB[i][7] =
dint(ZFFBU);
1186 BU_TAB[i][8] =
dint(AFFBU);
1187 BU_TAB[i][11] = NbLamH;
1188
1189
1190
1191
1192 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6],
1193 &VXOUT, &VYOUT, &VZOUT);
1194 BU_TAB[i][4] = VXOUT;
1195 BU_TAB[i][5] = VYOUT;
1196 BU_TAB[i][6] = VZOUT;
1197
1200 opt->optimfallowed = 0;
1201 fiss->ifis = 0;
1202
1203 G4double zffimf, affimf, zdummy1, adummy1, tkedummy1, jprf2, vx2ev_imf, vy2ev_imf,
1204 vz2ev_imf;
1205
1206 evapora(ZIMFBU, AIMFBU, &EEIMF, JPRFLIGHT, &zffimf, &affimf, &mtota, &vz2ev_imf,
1207 &vx2ev_imf, &vy2ev_imf, &FFBU2, &FIMFBU2, &zdummy1, &adummy1, &tkedummy1, &jprf2,
1208 &inttype, &inum, EV_TEMP, &IEV_TAB_TEMP, &NbLamimf);
1209
1210 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++) {
1211 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1212 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1213 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1214
1215
1216
1217
1218
1219
1220 lorentz_boost(BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], EV_TEMP[IJ][2], EV_TEMP[IJ][3],
1221 EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1222 lorentz_boost(VX2_IMF, VY2_IMF, VZ2_IMF, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1223 &VZ2OUT);
1224 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
1225 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
1226 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
1227 }
1228 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1229
1230 BU_TAB[IMULTBU + ILOOPBU][0] = BU_TAB[i][0];
1231 BU_TAB[IMULTBU + ILOOPBU][1] = BU_TAB[i][1];
1232 BU_TAB[IMULTBU + ILOOPBU][2] = BU_TAB[i][2];
1233 BU_TAB[IMULTBU + ILOOPBU][3] = BU_TAB[i][3];
1234 BU_TAB[IMULTBU + ILOOPBU][7] =
dint(zffimf);
1235 BU_TAB[IMULTBU + ILOOPBU][8] =
dint(affimf);
1236 BU_TAB[IMULTBU + ILOOPBU][11] = NbLamimf;
1237
1238 lorentz_boost(VX2_IMF, VY2_IMF, VZ2_IMF, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VXOUT,
1239 &VYOUT, &VZOUT);
1240 lorentz_boost(vx2ev_imf, vy2ev_imf, vz2ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1241 &VZ2OUT);
1242 BU_TAB[IMULTBU + ILOOPBU][4] = VX2OUT;
1243 BU_TAB[IMULTBU + ILOOPBU][5] = VY2OUT;
1244 BU_TAB[IMULTBU + ILOOPBU][6] = VZ2OUT;
1245 ILOOPBU = ILOOPBU + 1;
1246 }
1247 }
1248 else {
1249
1250
1251
1252
1253 BU_TAB[i][7] = BU_TAB[i][0];
1254 BU_TAB[i][8] = BU_TAB[i][1];
1255
1256
1257
1258 BU_TAB[i][11] = Nblamb[i];
1259 }
1260 }
1261
1262 IMULTBU = IMULTBU + ILOOPBU;
1263
1264
1265
1266 INEWLOOP = 0;
1267 ABU_SUM = 0.0;
1268 ZBU_SUM = 0.0;
1269
1270 for (
G4int i = 0; i < IMULTBU; i++) {
1271 ABU_SUM = ABU_SUM + BU_TAB[i][8];
1272 ZBU_SUM = ZBU_SUM + BU_TAB[i][7];
1274 BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VP1X, &VP1Y, &VP1Z, BU_TAB_TEMP1,
1275 &ILOOP);
1276
1277
1278
1279
1280
1281
1282 if (IOUNSTABLE > 0) {
1283
1284 ABU_SUM = ABU_SUM +
G4double(afpnew) - BU_TAB[i][8];
1285 ZBU_SUM = ZBU_SUM +
G4double(zfpnew) - BU_TAB[i][7];
1288 BU_TAB[i][4] = VP1X;
1289 BU_TAB[i][5] = VP1Y;
1290 BU_TAB[i][6] = VP1Z;
1291
1292
1293 for (
G4int IJ = 0; IJ < ILOOP; IJ++) {
1294 BU_TAB[IMULTBU + INEWLOOP + IJ][7] = BU_TAB_TEMP1[IJ][0];
1295 BU_TAB[IMULTBU + INEWLOOP + IJ][8] = BU_TAB_TEMP1[IJ][1];
1296 BU_TAB[IMULTBU + INEWLOOP + IJ][4] = BU_TAB_TEMP1[IJ][2];
1297 BU_TAB[IMULTBU + INEWLOOP + IJ][5] = BU_TAB_TEMP1[IJ][3];
1298 BU_TAB[IMULTBU + INEWLOOP + IJ][6] = BU_TAB_TEMP1[IJ][4];
1299 BU_TAB[IMULTBU + INEWLOOP + IJ][2] = 0.0;
1300 BU_TAB[IMULTBU + INEWLOOP + IJ][3] = 0.0;
1301 BU_TAB[IMULTBU + INEWLOOP + IJ][0] = BU_TAB[i][0];
1302 BU_TAB[IMULTBU + INEWLOOP + IJ][1] = BU_TAB[i][1];
1303 BU_TAB[IMULTBU + INEWLOOP + IJ][11] = BU_TAB[i][11];
1304 ABU_SUM = ABU_SUM + BU_TAB[IMULTBU + INEWLOOP + IJ][8];
1305 ZBU_SUM = ZBU_SUM + BU_TAB[IMULTBU + INEWLOOP + IJ][7];
1306 }
1307
1308 INEWLOOP = INEWLOOP + ILOOP;
1309 }
1310 }
1311
1312
1313 IMULTBU = IMULTBU + INEWLOOP;
1314
1315
1316 lorentz_boost(VX_incl, VY_incl, VZ_incl, VX_PREF, VY_PREF, VZ_PREF, &VXOUT, &VYOUT, &VZOUT);
1317 VX_PREF = VXOUT;
1318 VY_PREF = VYOUT;
1319 VZ_PREF = VZOUT;
1320
1321 for (
G4int i = 0; i < IMULTBU; i++) {
1322 lorentz_boost(VX_incl, VY_incl, VZ_incl, BU_TAB[i][4], BU_TAB[i][5], BU_TAB[i][6], &VXOUT,
1323 &VYOUT, &VZOUT);
1324 BU_TAB[i][4] = VXOUT;
1325 BU_TAB[i][5] = VYOUT;
1326 BU_TAB[i][6] = VZOUT;
1327 }
1328 for (
G4int i = 0; i < IEV_TAB; i++) {
1329 lorentz_boost(VX_incl, VY_incl, VZ_incl, EV_TAB[i][2], EV_TAB[i][3], EV_TAB[i][4], &VXOUT,
1330 &VYOUT, &VZOUT);
1331 EV_TAB[i][2] = VXOUT;
1332 EV_TAB[i][3] = VYOUT;
1333 EV_TAB[i][4] = VZOUT;
1334 }
1335 if (IMULTBU > 200) std::cout << "IMULTBU>200 " << IMULTBU << std::endl;
1336 }
1337
1338mult7777:
1339
1340
1343
1344 if (IMULTIFR == 0) {
1345
1346
1347 VX_PREF = VX_incl;
1348 VY_PREF = VY_incl;
1349 VZ_PREF = VZ_incl;
1350 }
1351
1352 if (IMULTIFR == 1) {
1353 NbLam0 = NbLamprf;
1354 }
1355
1356
1357
1358 opt->optimfallowed = 1;
1359 fiss->ifis = 1;
1360 fimf = 0;
1361 ff = 0;
1362
1363
1364
1365 if (zprfp <= 2 && zprfp < aprfp) {
1366 zf = zprf;
1367 af = aprf;
1368 ee = 0.0;
1369 ff = 0;
1370 fimf = 0;
1371 ftype = 0;
1372 aimf = 0.0;
1373 zimf = 0.0;
1374 tkeimf = 0.0;
1375 vx_eva = 0.0;
1376 vy_eva = 0.0;
1377 vz_eva = 0.0;
1378 jprf0 = jprf;
1379 goto a1972;
1380 }
1381
1382
1383 if (zprfp <= 2 && zprfp == aprfp) {
1384 unstable_nuclei(aprfp, zprfp, &afpnew, &zfpnew, IOUNSTABLE, VX_PREF, VY_PREF, VZ_PREF, &VP1X,
1385 &VP1Y, &VP1Z, EV_TAB_TEMP, &ILOOP);
1388 VX_PREF = VP1X;
1389 VY_PREF = VP1Y;
1390 VZ_PREF = VP1Z;
1391 for (
G4int I = 0; I < ILOOP; I++) {
1392 for (
G4int IJ = 0; IJ < 6; IJ++)
1393 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1394 }
1395 IEV_TAB = IEV_TAB + ILOOP;
1396 ee = 0.0;
1397 ff = 0;
1398 fimf = 0;
1399 ftype = 0;
1400 aimf = 0.0;
1401 zimf = 0.0;
1402 tkeimf = 0.0;
1403 vx_eva = 0.0;
1404 vy_eva = 0.0;
1405 vz_eva = 0.0;
1406 jprf0 = jprf;
1407 goto a1972;
1408 }
1409
1410
1411 if (zprfp == aprfp) {
1412 unstable_nuclei(aprfp, zprfp, &afpnew, &zfpnew, IOUNSTABLE, VX_PREF, VY_PREF, VZ_PREF, &VP1X,
1413 &VP1Y, &VP1Z, EV_TAB_TEMP, &ILOOP);
1416 VX_PREF = VP1X;
1417 VY_PREF = VP1Y;
1418 VZ_PREF = VP1Z;
1419 for (
G4int I = 0; I < ILOOP; I++) {
1420 for (
G4int IJ = 0; IJ < 6; IJ++)
1421 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1422 }
1423 IEV_TAB = IEV_TAB + ILOOP;
1424 ee = 0.0;
1425 ff = 0;
1426 fimf = 0;
1427 ftype = 0;
1428 aimf = 0.0;
1429 zimf = 0.0;
1430 tkeimf = 0.0;
1431 vx_eva = 0.0;
1432 vy_eva = 0.0;
1433 vz_eva = 0.0;
1434 jprf0 = jprf;
1435 goto a1972;
1436 }
1437
1438 evapora(zprf, aprf, &ee, jprf, &zf, &af, &mtota, &vz_eva, &vx_eva, &vy_eva, &ff, &fimf, &zimf,
1439 &aimf, &tkeimf, &jprf0, &inttype, &inum, EV_TEMP, &IEV_TAB_TEMP, &NbLam0);
1440
1441 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++) {
1442 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1443 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1444 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1445
1446
1447
1448
1449
1450 lorentz_boost(VX_PREF, VY_PREF, VZ_PREF, EV_TEMP[IJ][2], EV_TEMP[IJ][3], EV_TEMP[IJ][4], &VXOUT,
1451 &VYOUT, &VZOUT);
1452 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1453 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1454 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1455 }
1456 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1457
1458a1972:
1459
1460
1461 lorentz_boost(VX_PREF, VY_PREF, VZ_PREF, vx_eva, vy_eva, vz_eva, &VXOUT, &VYOUT, &VZOUT);
1462 V_CM[0] = VXOUT;
1463 V_CM[1] = VYOUT;
1464 V_CM[2] = VZOUT;
1465
1466 if (ff == 0 && fimf == 0) {
1467
1468 ftype = 0;
1471 SFP1 = NbLam0;
1472 AFPIMF = 0;
1473 ZFPIMF = 0;
1474 SFPIMF = 0;
1475 ZFP2 = 0;
1476 AFP2 = 0;
1477 SFP2 = 0;
1478 VFP1_CM[0] = V_CM[0];
1479 VFP1_CM[1] = V_CM[1];
1480 VFP1_CM[2] = V_CM[2];
1481 for (
G4int j = 0; j < 3; j++) {
1482 VIMF_CM[j] = 0.0;
1483 VFP2_CM[j] = 0.0;
1484 }
1485 }
1486
1487 if (ff == 1 && fimf == 0) ftype = 1;
1488 if (ff == 0 && fimf == 1)
1489 ftype = 2;
1490
1491
1492
1493
1494
1495
1496
1497 if (ftype == 1) {
1498 varntp->kfis = 1;
1499 if (NbLam0 > 0) varntp->kfis = 20;
1500
1501
1502 G4int IEV_TAB_FIS = 0, imode = 0;
1503
1504 G4double vx1_fission = 0., vy1_fission = 0., vz1_fission = 0.;
1505 G4double vx2_fission = 0., vy2_fission = 0., vz2_fission = 0.;
1506 G4double vx_eva_sc = 0., vy_eva_sc = 0., vz_eva_sc = 0.;
1507
1508 fission(af, zf, ee, jprf0, &vx1_fission, &vy1_fission, &vz1_fission, &vx2_fission, &vy2_fission,
1509 &vz2_fission, &ZFP1, &AFP1, &SFP1, &ZFP2, &AFP2, &SFP2, &imode, &vx_eva_sc, &vy_eva_sc,
1510 &vz_eva_sc, EV_TEMP, &IEV_TAB_FIS, &NbLam0);
1511
1512 for (
G4int IJ = 0; IJ < IEV_TAB_FIS; IJ++) {
1513 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1514 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1515 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1516
1517
1518
1519
1520
1521 lorentz_boost(V_CM[0], V_CM[1], V_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3], EV_TEMP[IJ][4],
1522 &VXOUT, &VYOUT, &VZOUT);
1523 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1524 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1525 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1526 }
1527 IEV_TAB = IEV_TAB + IEV_TAB_FIS;
1528
1529
1530
1531
1532 AFPIMF = 0;
1533 ZFPIMF = 0;
1534 SFPIMF = 0;
1535
1536
1537
1538
1539
1540
1541
1542
1543 lorentz_boost(vx1_fission, vy1_fission, vz1_fission, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT,
1544 &VZOUT);
1545 lorentz_boost(vx_eva_sc, vy_eva_sc, vz_eva_sc, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1546 VFP1_CM[0] = VX2OUT;
1547 VFP1_CM[1] = VY2OUT;
1548 VFP1_CM[2] = VZ2OUT;
1549
1550
1551
1552
1553
1554
1555 lorentz_boost(vx2_fission, vy2_fission, vz2_fission, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT,
1556 &VZOUT);
1557 lorentz_boost(vx_eva_sc, vy_eva_sc, vz_eva_sc, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1558 VFP2_CM[0] = VX2OUT;
1559 VFP2_CM[1] = VY2OUT;
1560 VFP2_CM[2] = VZ2OUT;
1561
1562
1563
1564
1565 }
1566 else if (ftype == 2) {
1567
1568
1571 opt->optimfallowed = 1;
1572 fiss->ifis = 1;
1573
1576 G4double pbH = (af - zf) / (af - zf + aimf - zimf);
1577
1578 for (
G4int i = 0; i < NbLam0; i++) {
1580 NbLamH++;
1581 }
1582 else {
1583 NbLamimf++;
1584 }
1585 }
1586
1587
1588 G4double EkinR1 = tkeimf * aimf / (af + aimf);
1589 G4double EkinR2 = tkeimf * af / (af + aimf);
1590 G4double V1 = std::sqrt(EkinR1 / af) * 1.3887;
1591 G4double V2 = std::sqrt(EkinR2 / aimf) * 1.3887;
1593 G4double VPERP1 = std::sqrt(
V1 *
V1 - VZ1_IMF * VZ1_IMF);
1595 G4double VX1_IMF = VPERP1 * std::sin(ALPHA1);
1596 G4double VY1_IMF = VPERP1 * std::cos(ALPHA1);
1600
1601 G4double EEIMFP = ee * af / (af + aimf);
1602 G4double EEIMF = ee * aimf / (af + aimf);
1603
1604
1606 0.40 * 931.490 * 1.160 * 1.160 * (std::pow(aimf, 5.0 / 3.0) + std::pow(af, 5.0 / 3.0))
1607 + 931.490 * 1.160 * 1.160 * aimf * af / (aimf + af)
1608 * (std::pow(aimf, 1. / 3.) + std::pow(af, 1. / 3.))
1609 * (std::pow(aimf, 1. / 3.) + std::pow(af, 1. / 3.));
1610
1611 G4double JPRFHEAVY = jprf0 * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(af, 5.0 / 3.0) / IINERTTOT;
1612 G4double JPRFLIGHT = jprf0 * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(aimf, 5.0 / 3.0) / IINERTTOT;
1613 if (af < 2.0) std::cout << "RN117-4,AF,ZF,EE,JPRFheavy" << std::endl;
1614
1615 G4double vx1ev_imf = 0., vy1ev_imf = 0., vz1ev_imf = 0., zdummy = 0., adummy = 0.,
1616 tkedummy = 0., jprf1 = 0.;
1617
1618 evapora(zf, af, &EEIMFP, JPRFHEAVY, &zff, &aff, &mtota, &vz1ev_imf, &vx1ev_imf, &vy1ev_imf,
1619 &FF11, &FIMF11, &zdummy, &adummy, &tkedummy, &jprf1, &inttype, &inum, EV_TEMP,
1620 &IEV_TAB_TEMP, &NbLamH);
1621
1622 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++) {
1623 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1624 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1625 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1626
1627
1628
1629
1630
1631 lorentz_boost(V_CM[0], V_CM[1], V_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3], EV_TEMP[IJ][4],
1632 &VXOUT, &VYOUT, &VZOUT);
1633 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1634 &VZ2OUT);
1635 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
1636 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
1637 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
1638 }
1639 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1640
1641
1644 opt->optimfallowed = 0;
1645 fiss->ifis = 0;
1646
1647
1648 G4double zffimf, affimf, zdummy1 = 0., adummy1 = 0., tkedummy1 = 0., jprf2, vx2ev_imf,
1649 vy2ev_imf, vz2ev_imf;
1650
1651 evapora(zimf, aimf, &EEIMF, JPRFLIGHT, &zffimf, &affimf, &mtota, &vz2ev_imf, &vx2ev_imf,
1652 &vy2ev_imf, &FF22, &FIMF22, &zdummy1, &adummy1, &tkedummy1, &jprf2, &inttype, &inum,
1653 EV_TEMP, &IEV_TAB_TEMP, &NbLamimf);
1654
1655 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++) {
1656 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1657 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1658 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1659
1660
1661
1662
1663
1664 lorentz_boost(V_CM[0], V_CM[1], V_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3], EV_TEMP[IJ][4],
1665 &VXOUT, &VYOUT, &VZOUT);
1666 lorentz_boost(VX2_IMF, VY2_IMF, VZ2_IMF, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1667 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
1668 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
1669 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
1670 }
1671 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1672
1673
1676 SFPIMF = NbLamimf;
1677
1678
1679
1680
1681
1682
1683
1684 lorentz_boost(VX2_IMF, VY2_IMF, VZ2_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
1685 lorentz_boost(vx2ev_imf, vy2ev_imf, vz2ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1686 VIMF_CM[0] = VX2OUT;
1687 VIMF_CM[1] = VY2OUT;
1688 VIMF_CM[2] = VZ2OUT;
1689
1690
1691
1692
1693 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
1694 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT, &VZ2OUT);
1695 VFP1_CM[0] = VX2OUT;
1696 VFP1_CM[1] = VY2OUT;
1697 VFP1_CM[2] = VZ2OUT;
1698
1699 if (FF11 == 0 && FIMF11 == 0) {
1700
1703 SFP1 = NbLamH;
1704 ZFP2 = 0;
1705 AFP2 = 0;
1706 SFP2 = 0;
1707 ftype = 2;
1710 SFPIMF = NbLamimf;
1711 for (
G4int I = 0; I < 3; I++)
1712 VFP2_CM[I] = 0.0;
1713 }
1714 else if (FF11 == 1 && FIMF11 == 0) {
1715
1716 varntp->kfis = 1;
1717 if (NbLam0 > 0) varntp->kfis = 20;
1718
1719 opt->optimfallowed = 0;
1720 fiss->ifis = 0;
1721
1722 zf = zff;
1723 af = aff;
1724 ee = EEIMFP;
1725
1726 ftype = 21;
1727
1728 G4int IEV_TAB_FIS = 0, imode = 0;
1729
1730 G4double vx1_fission = 0., vy1_fission = 0., vz1_fission = 0.;
1731 G4double vx2_fission = 0., vy2_fission = 0., vz2_fission = 0.;
1732 G4double vx_eva_sc = 0., vy_eva_sc = 0., vz_eva_sc = 0.;
1733
1734 fission(af, zf, ee, jprf1, &vx1_fission, &vy1_fission, &vz1_fission, &vx2_fission,
1735 &vy2_fission, &vz2_fission, &ZFP1, &AFP1, &SFP1, &ZFP2, &AFP2, &SFP2, &imode,
1736 &vx_eva_sc, &vy_eva_sc, &vz_eva_sc, EV_TEMP, &IEV_TAB_FIS, &NbLamH);
1737
1738 for (
G4int IJ = 0; IJ < IEV_TAB_FIS; IJ++) {
1739 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1740 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1741 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1742
1743
1744
1745
1746
1747 lorentz_boost(VFP1_CM[0], VFP1_CM[1], VFP1_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3],
1748 EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1749 EV_TAB[IJ + IEV_TAB][2] = VXOUT;
1750 EV_TAB[IJ + IEV_TAB][3] = VYOUT;
1751 EV_TAB[IJ + IEV_TAB][4] = VZOUT;
1752 }
1753 IEV_TAB = IEV_TAB + IEV_TAB_FIS;
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
1766 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1767 &VZ2OUT);
1768 lorentz_boost(vx1_fission, vy1_fission, vz1_fission, VX2OUT, VY2OUT, VZ2OUT, &VXOUT, &VYOUT,
1769 &VZOUT);
1770 lorentz_boost(vx_eva_sc, vy_eva_sc, vz_eva_sc, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1771 &VZ2OUT);
1772 VFP1_CM[0] = VX2OUT;
1773 VFP1_CM[1] = VY2OUT;
1774 VFP1_CM[2] = VZ2OUT;
1775
1776
1777
1778
1779
1780
1781
1782
1783 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
1784 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1785 &VZ2OUT);
1786 lorentz_boost(vx2_fission, vy2_fission, vz2_fission, VX2OUT, VY2OUT, VZ2OUT, &VXOUT, &VYOUT,
1787 &VZOUT);
1788 lorentz_boost(vx_eva_sc, vy_eva_sc, vz_eva_sc, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1789 &VZ2OUT);
1790 VFP2_CM[0] = VX2OUT;
1791 VFP2_CM[1] = VY2OUT;
1792 VFP2_CM[2] = VZ2OUT;
1793 }
1794 else if (FF11 == 0 && FIMF11 == 1) {
1795
1796
1797 opt->optimfallowed = 0;
1798 fiss->ifis = 0;
1799
1800 zf = zff;
1801 af = aff;
1802 ee = EEIMFP;
1803 aimf = adummy;
1804 zimf = zdummy;
1805 tkeimf = tkedummy;
1806 FF11 = 0;
1807 FIMF11 = 0;
1808 ftype = 22;
1809
1811 G4int NbLamimf1 = 0;
1812 G4double pbH1 = (af - zf) / (af - zf + aimf - zimf);
1813 for (
G4int i = 0; i < NbLamH; i++) {
1815 NbLamH1++;
1816 }
1817 else {
1818 NbLamimf1++;
1819 }
1820 }
1821
1822
1823 EkinR1 = tkeimf * aimf / (af + aimf);
1824 EkinR2 = tkeimf * af / (af + aimf);
1825 V1 = std::sqrt(EkinR1 / af) * 1.3887;
1826 V2 = std::sqrt(EkinR2 / aimf) * 1.3887;
1828 VPERP1 = std::sqrt(
V1 *
V1 - VZ1_IMFS * VZ1_IMFS);
1830 G4double VX1_IMFS = VPERP1 * std::sin(ALPHA1);
1831 G4double VY1_IMFS = VPERP1 * std::cos(ALPHA1);
1835
1836 EEIMFP = ee * af / (af + aimf);
1837 EEIMF = ee * aimf / (af + aimf);
1838
1839
1840 IINERTTOT =
1841 0.40 * 931.490 * 1.160 * 1.160 * (std::pow(aimf, 5.0 / 3.0) + std::pow(af, 5.0 / 3.0))
1842 + 931.490 * 1.160 * 1.160 * aimf * af / (aimf + af)
1843 * (std::pow(aimf, 1. / 3.) + std::pow(af, 1. / 3.))
1844 * (std::pow(aimf, 1. / 3.) + std::pow(af, 1. / 3.));
1845
1846 JPRFHEAVY = jprf1 * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(af, 5.0 / 3.0) / IINERTTOT;
1847 JPRFLIGHT = jprf1 * 0.4 * 931.49 * 1.16 * 1.16 * std::pow(aimf, 5.0 / 3.0) / IINERTTOT;
1848
1849 G4double zffs = 0., affs = 0., vx1ev_imfs = 0., vy1ev_imfs = 0., vz1ev_imfs = 0., jprf3 = 0.;
1850
1851 evapora(zf, af, &EEIMFP, JPRFHEAVY, &zffs, &affs, &mtota, &vz1ev_imfs, &vx1ev_imfs,
1852 &vy1ev_imfs, &FF11, &FIMF11, &zdummy, &adummy, &tkedummy, &jprf3, &inttype, &inum,
1853 EV_TEMP, &IEV_TAB_TEMP, &NbLamH1);
1854
1855 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++) {
1856 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1857 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1858 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1859
1860
1861
1862
1863
1864 lorentz_boost(VFP1_CM[0], VFP1_CM[1], VFP1_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3],
1865 EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1866 lorentz_boost(vx1ev_imfs, vy1ev_imfs, vz1ev_imfs, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1867 &VZ2OUT);
1868 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
1869 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
1870 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
1871 }
1872 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1873
1874
1875 opt->optimfallowed = 0;
1876 fiss->ifis = 0;
1877
1878 FF22 = 0;
1879 FIMF22 = 0;
1880
1881 G4double zffimfs = 0., affimfs = 0., vx2ev_imfs = 0., vy2ev_imfs = 0., vz2ev_imfs = 0.,
1882 jprf4 = 0.;
1883
1884 evapora(zimf, aimf, &EEIMF, JPRFLIGHT, &zffimfs, &affimfs, &mtota, &vz2ev_imfs, &vx2ev_imfs,
1885 &vy2ev_imfs, &FF22, &FIMF22, &zdummy1, &adummy1, &tkedummy1, &jprf4, &inttype, &inum,
1886 EV_TEMP, &IEV_TAB_TEMP, &NbLamimf1);
1887
1888 for (
G4int IJ = 0; IJ < IEV_TAB_TEMP; IJ++) {
1889 EV_TAB[IJ + IEV_TAB][0] = EV_TEMP[IJ][0];
1890 EV_TAB[IJ + IEV_TAB][1] = EV_TEMP[IJ][1];
1891 EV_TAB[IJ + IEV_TAB][5] = EV_TEMP[IJ][5];
1892
1893
1894
1895
1896
1897 lorentz_boost(VFP1_CM[0], VFP1_CM[1], VFP1_CM[2], EV_TEMP[IJ][2], EV_TEMP[IJ][3],
1898 EV_TEMP[IJ][4], &VXOUT, &VYOUT, &VZOUT);
1899 lorentz_boost(vx2ev_imfs, vy2ev_imfs, vz2ev_imfs, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1900 &VZ2OUT);
1901 EV_TAB[IJ + IEV_TAB][2] = VX2OUT;
1902 EV_TAB[IJ + IEV_TAB][3] = VY2OUT;
1903 EV_TAB[IJ + IEV_TAB][4] = VZ2OUT;
1904 }
1905 IEV_TAB = IEV_TAB + IEV_TAB_TEMP;
1906
1909 SFP1 = NbLamH1;
1912 SFP2 = NbLamimf1;
1913
1914
1915
1916
1917
1918
1919 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
1920 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1921 &VZ2OUT);
1922 lorentz_boost(VX1_IMFS, VY1_IMFS, VZ1_IMFS, VX2OUT, VY2OUT, VZ2OUT, &VXOUT, &VYOUT, &VZOUT);
1923 lorentz_boost(vx1ev_imfs, vy1ev_imfs, vz1ev_imfs, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1924 &VZ2OUT);
1925 VFP1_CM[0] = VX2OUT;
1926 VFP1_CM[1] = VY2OUT;
1927 VFP1_CM[2] = VZ2OUT;
1928
1929
1930
1931
1932
1933
1934 lorentz_boost(VX1_IMF, VY1_IMF, VZ1_IMF, V_CM[0], V_CM[1], V_CM[2], &VXOUT, &VYOUT, &VZOUT);
1935 lorentz_boost(vx1ev_imf, vy1ev_imf, vz1ev_imf, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1936 &VZ2OUT);
1937 lorentz_boost(VX2_IMFS, VY2_IMFS, VZ2_IMFS, VX2OUT, VY2OUT, VZ2OUT, &VXOUT, &VYOUT, &VZOUT);
1938 lorentz_boost(vx2ev_imfs, vy2ev_imfs, vz2ev_imfs, VXOUT, VYOUT, VZOUT, &VX2OUT, &VY2OUT,
1939 &VZ2OUT);
1940 VFP2_CM[0] = VX2OUT;
1941 VFP2_CM[1] = VY2OUT;
1942 VFP2_CM[2] = VZ2OUT;
1943 }
1944 }
1945
1946
1947 if (ftype != 1 && ftype != 21) {
1948
1949 IOUNSTABLE = 0;
1950
1951 unstable_nuclei(AFP1, ZFP1, &afpnew, &zfpnew, IOUNSTABLE, VFP1_CM[0], VFP1_CM[1], VFP1_CM[2],
1952 &VP1X, &VP1Y, &VP1Z, EV_TAB_TEMP, &ILOOP);
1953
1954 if (IOUNSTABLE == 1) {
1955 AFP1 = afpnew;
1956 ZFP1 = zfpnew;
1957 VFP1_CM[0] = VP1X;
1958 VFP1_CM[1] = VP1Y;
1959 VFP1_CM[2] = VP1Z;
1960 for (
G4int I = 0; I < ILOOP; I++) {
1961 for (
G4int IJ = 0; IJ < 5; IJ++)
1962 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1963 }
1964 IEV_TAB = IEV_TAB + ILOOP;
1965 }
1966
1967 if (ftype > 1) {
1968 IOUNSTABLE = 0;
1969
1970 unstable_nuclei(AFPIMF, ZFPIMF, &afpnew, &zfpnew, IOUNSTABLE, VIMF_CM[0], VIMF_CM[1],
1971 VIMF_CM[2], &VP1X, &VP1Y, &VP1Z, EV_TAB_TEMP, &ILOOP);
1972
1973 if (IOUNSTABLE == 1) {
1974 AFPIMF = afpnew;
1975 ZFPIMF = zfpnew;
1976 VIMF_CM[0] = VP1X;
1977 VIMF_CM[1] = VP1Y;
1978 VIMF_CM[2] = VP1Z;
1979 for (
G4int I = 0; I < ILOOP; I++) {
1980 for (
G4int IJ = 0; IJ < 5; IJ++)
1981 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
1982 }
1983 IEV_TAB = IEV_TAB + ILOOP;
1984 }
1985
1986 if (ftype > 2) {
1987 IOUNSTABLE = 0;
1988
1989 unstable_nuclei(AFP2, ZFP2, &afpnew, &zfpnew, IOUNSTABLE, VFP2_CM[0], VFP2_CM[1],
1990 VFP2_CM[2], &VP1X, &VP1Y, &VP1Z, EV_TAB_TEMP, &ILOOP);
1991
1992 if (IOUNSTABLE == 1) {
1993 AFP2 = afpnew;
1994 ZFP2 = zfpnew;
1995 VFP2_CM[0] = VP1X;
1996 VFP2_CM[1] = VP1Y;
1997 VFP2_CM[2] = VP1Z;
1998 for (
G4int I = 0; I < ILOOP; I++) {
1999 for (
G4int IJ = 0; IJ < 5; IJ++)
2000 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2001 }
2002 IEV_TAB = IEV_TAB + ILOOP;
2003 }
2004 }
2005 }
2006 }
2007
2008
2009 if (ftype == 1 || ftype == 21) {
2010
2011 IOUNSTABLE = 0;
2012
2013 unstable_nuclei(AFP1, ZFP1, &afpnew, &zfpnew, IOUNSTABLE, VFP1_CM[0], VFP1_CM[1], VFP1_CM[2],
2014 &VP1X, &VP1Y, &VP1Z, EV_TAB_TEMP, &ILOOP);
2015
2016 if (IOUNSTABLE == 1) {
2017 AFP1 = afpnew;
2018 ZFP1 = zfpnew;
2019 VFP1_CM[0] = VP1X;
2020 VFP1_CM[1] = VP1Y;
2021 VFP1_CM[2] = VP1Z;
2022 for (
G4int I = 0; I < ILOOP; I++) {
2023 for (
G4int IJ = 0; IJ < 5; IJ++)
2024 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2025 }
2026 IEV_TAB = IEV_TAB + ILOOP;
2027 }
2028
2029 IOUNSTABLE = 0;
2030
2031 unstable_nuclei(AFP2, ZFP2, &afpnew, &zfpnew, IOUNSTABLE, VFP2_CM[0], VFP2_CM[1], VFP2_CM[2],
2032 &VP1X, &VP1Y, &VP1Z, EV_TAB_TEMP, &ILOOP);
2033
2034 if (IOUNSTABLE == 1) {
2035 AFP2 = afpnew;
2036 ZFP2 = zfpnew;
2037 VFP2_CM[0] = VP1X;
2038 VFP2_CM[1] = VP1Y;
2039 VFP2_CM[2] = VP1Z;
2040 for (
G4int I = 0; I < ILOOP; I++) {
2041 for (
G4int IJ = 0; IJ < 5; IJ++)
2042 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2043 }
2044 IEV_TAB = IEV_TAB + ILOOP;
2045 }
2046
2047 if (ftype == 21) {
2048 IOUNSTABLE = 0;
2049
2050 unstable_nuclei(AFPIMF, ZFPIMF, &afpnew, &zfpnew, IOUNSTABLE, VIMF_CM[0], VIMF_CM[1],
2051 VIMF_CM[2], &VP1X, &VP1Y, &VP1Z, EV_TAB_TEMP, &ILOOP);
2052
2053 if (IOUNSTABLE == 1) {
2054 AFPIMF = afpnew;
2055 ZFPIMF = zfpnew;
2056 VIMF_CM[0] = VP1X;
2057 VIMF_CM[1] = VP1Y;
2058 VIMF_CM[2] = VP1Z;
2059 for (
G4int I = 0; I < ILOOP; I++) {
2060 for (
G4int IJ = 0; IJ < 5; IJ++)
2061 EV_TAB[I + IEV_TAB][IJ] = EV_TAB_TEMP[I][IJ];
2062 }
2063 IEV_TAB = IEV_TAB + ILOOP;
2064 }
2065 }
2066 }
2067
2068
2069 if ((ftype == 1 || ftype == 21) && (AFP2 <= 0 || AFP1 <= 0 || ZFP2 <= 0 || ZFP1 <= 0)) {
2070 std::cout << "ZFP1:" << ZFP1 << std::endl;
2071 std::cout << "AFP1:" << AFP1 << std::endl;
2072 std::cout << "ZFP2:" << ZFP2 << std::endl;
2073 std::cout << "AFP2:" << AFP2 << std::endl;
2074 }
2075
2076
2077 EV_TAB[IEV_TAB][0] = ZFP1;
2078 EV_TAB[IEV_TAB][1] = AFP1;
2079 EV_TAB[IEV_TAB][5] = SFP1;
2080 EV_TAB[IEV_TAB][2] = VFP1_CM[0];
2081 EV_TAB[IEV_TAB][3] = VFP1_CM[1];
2082 EV_TAB[IEV_TAB][4] = VFP1_CM[2];
2083 IEV_TAB = IEV_TAB + 1;
2084
2085 if (AFP2 > 0) {
2086 EV_TAB[IEV_TAB][0] = ZFP2;
2087 EV_TAB[IEV_TAB][1] = AFP2;
2088 EV_TAB[IEV_TAB][5] = SFP2;
2089 EV_TAB[IEV_TAB][2] = VFP2_CM[0];
2090 EV_TAB[IEV_TAB][3] = VFP2_CM[1];
2091 EV_TAB[IEV_TAB][4] = VFP2_CM[2];
2092 IEV_TAB = IEV_TAB + 1;
2093 }
2094
2095 if (AFPIMF > 0) {
2096 EV_TAB[IEV_TAB][0] = ZFPIMF;
2097 EV_TAB[IEV_TAB][1] = AFPIMF;
2098 EV_TAB[IEV_TAB][5] = SFPIMF;
2099 EV_TAB[IEV_TAB][2] = VIMF_CM[0];
2100 EV_TAB[IEV_TAB][3] = VIMF_CM[1];
2101 EV_TAB[IEV_TAB][4] = VIMF_CM[2];
2102 IEV_TAB = IEV_TAB + 1;
2103 }
2104
2106 return;
2107}
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.16.1