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