Changeset 1347 for trunk/source/processes/hadronic/models/incl/src
- Timestamp:
- Dec 22, 2010, 3:52:27 PM (14 years ago)
- Location:
- trunk/source/processes/hadronic/models/incl/src
- Files:
-
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/incl/src/G4Abla.cc
r1340 r1347 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4Abla.cc,v 1.2 2 2010/10/26 02:47:59 kaitanie Exp $26 // $Id: G4Abla.cc,v 1.27 2010/11/17 20:19:09 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 40 40 #include "G4AblaFissionSimfis18.hh" 41 41 #include "G4AblaFission.hh" 42 43 G4Abla::G4Abla()44 {45 verboseLevel = 0;46 ilast = 0;47 }48 49 G4Abla::G4Abla(G4Hazard *hazard, G4Volant *volant)50 {51 verboseLevel = 0;52 ilast = 0;53 volant = volant; // ABLA internal particle data54 volant->iv = 0;55 hazard = hazard; // Random seeds56 57 randomGenerator = new G4InclGeant4Random();58 //randomGenerator = new G4Ranecu();59 varntp = new G4VarNtp();60 pace = new G4Pace();61 ald = new G4Ald();62 eenuc = new G4Eenuc();63 ec2sub = new G4Ec2sub();64 ecld = new G4Ecld();65 fb = new G4Fb();66 fiss = new G4Fiss();67 opt = new G4Opt();68 }69 42 70 43 G4Abla::G4Abla(G4Hazard *aHazard, G4Volant *aVolant, G4VarNtp *aVarntp) … … 109 82 G4Abla::~G4Abla() 110 83 { 84 delete fissionModel; 111 85 delete randomGenerator; 112 86 delete pace; … … 483 457 // G4double *pleva_par, G4double *pxeva_par, G4double *pyeva_par, 484 458 // G4double *ff_par, G4int *inttype_par, G4int *inum_par); 485 G4double zf1 , af1, malpha1, ffpleva1, ffpxeva1, ffpyeva1;486 G4int ff1 , ftype1;459 G4double zf1 = 0.0, af1 = 0.0, malpha1 = 0.0, ffpleva1 = 0.0, ffpxeva1 = 0.0, ffpyeva1 = 0.0; 460 G4int ff1 = 0, ftype1 = 0; 487 461 evapora(zff1, aff1, epf1_out, 0.0, &zf1, &af1, &malpha1, &ffpleva1, 488 462 &ffpxeva1, &ffpyeva1, &ff1, &ftype1, &inum); … … 560 534 // G4double *pleva_par, G4double *pxeva_par, G4double *pyeva_par, 561 535 // G4double *ff_par, G4int *inttype_par, G4int *inum_par); 562 G4double zf2 , af2, malpha2, ffpleva2, ffpxeva2, ffpyeva2;563 G4int ff2 , ftype2;536 G4double zf2 = 0.0, af2 = 0.0, malpha2 = 0.0, ffpleva2 = 0.0, ffpxeva2 = 0.0, ffpyeva2 = 0.0; 537 G4int ff2 = 0, ftype2 = 0; 564 538 evapora(zff2,aff2,epf2_out,0.0,&zf2,&af2,&malpha2,&ffpleva2, 565 539 &ffpxeva2,&ffpyeva2,&ff2,&ftype2,&inum); … … 1608 1582 G4double *sbp_par, G4double *sba_par, G4double *ecn_par, 1609 1583 G4double *ecp_par,G4double *eca_par, G4double *bp_par, 1610 G4double *ba_par, G4int inttype, G4int inum, G4int itest)1584 G4double *ba_par, G4int, G4int inum, G4int itest) 1611 1585 { 1612 1586 G4int dummy0 = 0; … … 1758 1732 1759 1733 imaxwell = 1; 1760 inttype = 0;1761 1734 1762 1735 // limiting of excitation energy where fission occurs … … 2842 2815 i = idint(y/(2.0e-02)) + 1; 2843 2816 2844 if( i>= bsbkSize) {2817 if((i + 1) >= bsbkSize) { 2845 2818 if(verboseLevel > 2) { 2846 G4cout <<"G4Abla error: index i = " << i<< " is greater than array size permits." << G4endl;2819 G4cout <<"G4Abla error: index " << i + 1 << " is greater than array size permits." << G4endl; 2847 2820 } 2848 2821 bipolResult = 0.0; … … 3704 3677 // own class 3705 3678 3706 void G4Abla::standardRandom(G4double *rndm, G4long *seed) 3707 { 3708 (*seed) = (*seed); // Avoid warning during compilation. 3679 void G4Abla::standardRandom(G4double *rndm, G4long*) 3680 { 3709 3681 // Use Geant4 G4UniformRand 3710 3682 (*rndm) = randomGenerator->getRandom(); -
trunk/source/processes/hadronic/models/incl/src/G4AblaFission.cc
r968 r1347 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4AblaFission.cc,v 1. 3 2008/11/06 08:42:00 gcosmoExp $26 // $Id: G4AblaFission.cc,v 1.5 2010/11/17 20:19:09 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 36 36 G4AblaFission::G4AblaFission() 37 37 { 38 hazard = 0; 39 randomGenerator = 0; 38 40 } 39 41 … … 1086 1088 } 1087 1089 1088 void G4AblaFission::standardRandom(G4double *rndm, G4long *seed) 1089 { 1090 (*seed) = (*seed); // Avoid warning during compilation. 1090 void G4AblaFission::standardRandom(G4double *rndm, G4long*) 1091 { 1091 1092 // Use Geant4 G4UniformRand 1092 1093 (*rndm) = randomGenerator->getRandom(); -
trunk/source/processes/hadronic/models/incl/src/G4AblaFissionSimfis18.cc
r968 r1347 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4AblaFissionSimfis18.cc,v 1. 3 2008/11/06 08:42:00 gcosmoExp $26 // $Id: G4AblaFissionSimfis18.cc,v 1.5 2010/11/17 20:19:09 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 36 36 G4AblaFissionSimfis18::G4AblaFissionSimfis18() 37 37 { 38 hazard = 0; 39 randomGenerator = 0; 38 40 } 39 41 … … 1412 1414 } 1413 1415 1414 void G4AblaFissionSimfis18::standardRandom(G4double *rndm, G4long *seed) 1415 { 1416 (*seed) = (*seed); // Avoid warning during compilation. 1416 void G4AblaFissionSimfis18::standardRandom(G4double *rndm, G4long*) 1417 { 1417 1418 // Use Geant4 G4UniformRand 1418 1419 (*rndm) = randomGenerator->getRandom(); -
trunk/source/processes/hadronic/models/incl/src/G4Incl.cc
r1340 r1347 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4Incl.cc,v 1.3 3 2010/10/29 06:48:43gunter Exp $26 // $Id: G4Incl.cc,v 1.37 2010/12/15 07:41:31 gunter Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 57 57 randomGenerator = new G4InclGeant4Random(); 58 58 // randomGenerator = new G4Ranecu(); 59 clearState(); 59 60 } 60 61 … … 85 86 // randomGenerator = new G4Ranecu(); 86 87 randomGenerator = new G4InclGeant4Random(); 88 clearState(); 87 89 } 88 90 … … 149 151 150 152 theLogger = 0; 153 clearState(); 151 154 } 152 155 153 156 G4Incl::~G4Incl() 154 157 { 158 delete be; 159 delete ps; 160 delete fermi; 161 delete qvp; 162 155 163 delete randomGenerator; 156 164 delete light_gaus_nuc; … … 385 393 void G4Incl::processEventIncl(G4InclInput *input) 386 394 { 395 //G4cout <<"Starting event " << eventnumber << G4endl; 387 396 if(input == 0) { 388 397 G4cerr <<"G4Incl fatal error: NULL pointer passed as input!" << G4endl; … … 393 402 const G4double uma = 931.4942; 394 403 const G4double melec = 0.511; 395 396 G4double pcorem = 0.0; 397 G4double pxrem = 0.0; 398 G4double pyrem = 0.0; 399 G4double pzrem = 0.0; 404 const G4double fmp = 938.2796; 405 406 G4double pcorem = 0.0; 407 G4double pxrem = 0.0; 408 G4double pyrem = 0.0; 409 G4double pzrem = 0.0; 400 410 401 411 G4double ap = 0.0, zp = 0.0, mprojo = 0.0, pbeam = 0.0; 402 412 403 413 varntp->clear(); 404 if(calincl->bulletType() == 3) { // pi+ 414 415 if(calincl->bulletType() == -12) { 416 be->ia_be = std::abs(calincl->bulletType()); 417 be->iz_be = 6; 418 } else if(calincl->bulletType() == -666) { 419 be->iz_be = calincl->extendedProjectileZ(); 420 be->ia_be = calincl->extendedProjectileA(); 421 } 422 423 if(calincl->isExtendedProjectile() == false && calincl->bulletType() < -max_a_proj) { 424 // if(calincl->bulletType() < -max_a_proj) { 425 G4cout <<"max a of composite projectile is: " << max_a_proj << G4endl; 426 exit(0); 427 } 428 if(calincl->bulletType() < 0) { 429 // calincl->bulletType() = std::floor(calincl->bulletType() + 0.1); WTF??? 430 be->pms_be=100.; 431 G4int i_tabled=0; 432 if(be->iz_be == 3 && be->ia_be == 6) { 433 be->rms_be=2.56; 434 be->bind_be=32.0; 435 i_tabled=1; 436 } else if(be->iz_be == 3 && be->ia_be == 7) { // TODO: Check the values! 437 be->rms_be=2.56; 438 be->bind_be=32.0; 439 i_tabled=1; 440 } else if(be->iz_be == 3 && be->ia_be == 8) { 441 be->rms_be=2.40; 442 be->bind_be=39.25; 443 i_tabled=1; 444 } else if(be->iz_be == 4 && be->ia_be == 7) { 445 be->rms_be=2.51; 446 be->bind_be=58.17; 447 i_tabled=1; 448 } else if(be->iz_be == 4 && be->ia_be == 9) { 449 be->rms_be=2.51; 450 be->bind_be=58.17; 451 i_tabled=1; 452 } else if(be->iz_be == 4 && be->ia_be == 10) { 453 be->rms_be=2.45; 454 be->bind_be=64.75; 455 i_tabled=1; 456 } else if(be->iz_be == 5 && be->ia_be == 10) { 457 be->rms_be=2.45; 458 be->bind_be=64.75; 459 i_tabled=1; 460 } else if(be->iz_be == 5 && be->ia_be == 11) { 461 be->rms_be=2.40; 462 be->bind_be=76.21; 463 i_tabled=1; 464 } else if(be->iz_be == 6 && be->ia_be == 9) { // TODO: Check the values! 465 be->rms_be=2.44; 466 be->bind_be=92.17; 467 i_tabled=1; 468 } else if(be->iz_be == 6 && be->ia_be == 10) { // TODO: Check the values! 469 be->rms_be=2.44; 470 be->bind_be=92.17; 471 i_tabled=1; 472 } else if(be->iz_be == 6 && be->ia_be == 11) { // TODO: Check the values! 473 be->rms_be=2.44; 474 be->bind_be=92.17; 475 i_tabled=1; 476 } else if(be->iz_be == 6 && calincl->bulletType() == -12) { // Special Carbon case 477 G4cout <<"Carbon 12 (special) selected." << G4endl; 478 be->rms_be=2.44; 479 be->bind_be=92.17; 480 i_tabled=1; 481 } else if(be->iz_be == 6 && be->ia_be == 12) { 482 be->rms_be=2.44; 483 be->bind_be=92.17; 484 i_tabled=1; 485 } else if(be->iz_be == 7 && be->ia_be == 16) { 486 be->rms_be=2.73; 487 be->bind_be=127.62; 488 i_tabled=1; 489 } else { 490 G4cout <<"Warning: No rms and binding for projectile ion A = " << be->ia_be << " Z = " << be->iz_be << G4endl; 491 be->rms_be=2.44; 492 be->bind_be=92.17; 493 G4cout <<"Warning: Using probably bad values rms = " << be->rms_be << " binding = " << be->bind_be << G4endl; 494 i_tabled=1; 495 } 496 497 if(i_tabled == 0) { 498 G4cout <<"This heavy ion (a,z)= " << be->ia_be << " " << be->iz_be << " is not defined as beam in INCL" << G4endl; 499 exit(0); 500 } 501 502 // G4cout <<"z projectile, rms_r, rms_p (gaussian model)" << be->iz_be << " " << be->rms_be << " " << be->pms_be << G4endl; 503 // G4cout <<"binding energy (mev):" << be->bind_be << G4endl; 504 // G4cout <<"fermi-breakup dresner below a=" << calincl->f[11] << G4endl; 505 } 506 // G4cout <<"Target Mass and Charge: " << calincl->targetA() << " " << calincl->targetZ() << G4endl; 507 // calincl->f[10] = 0; // No clusters 508 509 if(calincl->bulletType() == -12) { // C12 special case 510 mprojo=fmp*std::abs(calincl->bulletType()) - be->bind_be; 511 pbeam=std::sqrt(calincl->bulletE()*(calincl->bulletE()+2.*mprojo)); 512 ap=std::abs(calincl->bulletType()); 513 zp=be->iz_be; 514 } else if(calincl->bulletType() == -666) { // Generic extended projectile 515 mprojo=fmp*be->ia_be - be->bind_be; 516 pbeam=std::sqrt(calincl->bulletE()*(calincl->bulletE()+2.*mprojo)); 517 ap=be->ia_be; 518 zp=be->iz_be; 519 } 520 // pi+ 521 if(calincl->bulletType() == 3) { 405 522 mprojo = 139.56995; 406 523 ap = 0.0; … … 408 525 } 409 526 410 if(calincl->bulletType() == 4) { // pi0 527 // pi0 528 if(calincl->bulletType() == 4) { 411 529 mprojo = 134.9764; 412 530 ap = 0.0; … … 414 532 } 415 533 416 if(calincl->bulletType() == 5) { // pi- 534 // pi- 535 if(calincl->bulletType() == 5) { 417 536 mprojo = 139.56995; 418 537 ap = 0.0; … … 420 539 } 421 540 422 // Coulomb en entree seulement pour les particules ci-dessous. 423 if(calincl->bulletType() == 1) { // proton 541 // coulomb en entree seulement pour les particules ci-dessous 542 543 // proton 544 if(calincl->bulletType() == 1) { 424 545 mprojo = 938.27231; 425 546 ap = 1.0; 426 547 zp = 1.0; 427 548 } 428 429 if(calincl->bulletType() == 2) { // neutron 549 550 // neutron 551 if(calincl->bulletType() == 2) { 430 552 mprojo = 939.56563; 431 553 ap = 1.0; 432 554 zp = 0.0; 433 555 } 434 435 if(calincl->bulletType() == 6) { // deuteron 556 557 // deuteron 558 if(calincl->bulletType() == 6) { 436 559 mprojo = 1875.61276; 437 560 ap = 2.0; 438 561 zp = 1.0; 439 562 } 440 441 if(calincl->bulletType() == 7) { // triton 563 564 // triton 565 if(calincl->bulletType() == 7) { 442 566 mprojo = 2808.95; 443 567 ap = 3.0; 444 568 zp = 1.0; 445 569 } 446 447 if(calincl->bulletType() == 8) { // He3 570 571 // He3 572 if(calincl->bulletType() == 8) { 448 573 mprojo = 2808.42; 449 574 ap = 3.0; … … 451 576 } 452 577 453 if(calincl->bulletType() == 9) { // Alpha 578 // Alpha 579 if(calincl->bulletType() == 9) { 454 580 mprojo = 3727.42; 455 581 ap = 4.0; … … 457 583 } 458 584 459 if(calincl->bulletType() == -12) { // Carbon 460 mprojo = 12.0 * 938.27231; 585 // Carbon 586 if(calincl->bulletType() == -12) { 587 mprojo = 6.0*938.27231 + 6.0*939.56563; 461 588 ap = 12.0; 462 589 zp = 6.0; 463 } else if(calincl->bulletType() == -666) { // Generic extended ion projectile464 mprojo = calincl->extendedProjectileA() * 938.27231;465 ap = calincl->extendedProjectileA();466 zp = calincl->extendedProjectileZ();467 590 } 468 591 … … 482 605 G4double bimpac = 0.0; 483 606 G4int jrem = 0; 607 G4double xjrem = 0.0, yjrem = 0.0, zjrem = 0.0; 484 608 G4double alrem = 0.0; 485 609 486 /** 487 * Coulomb barrier treatment. 488 */ 610 // Coulomb barrier 611 489 612 G4double probaTrans = 0.0; 490 613 G4double rndm = 0.0; 491 614 492 615 if((calincl->bulletType() == 1) || (calincl->bulletType() >= 6)) { 616 // probaTrans = coulombTransm(calincl->bulletE(),apro,zpro,calincl->targetA(),calincl->targetZ()); 493 617 probaTrans = coulombTransm(calincl->bulletE(),ap,zp,calincl->targetA(),calincl->targetZ()); 494 618 standardRandom(&rndm, &(hazard->ial)); … … 499 623 } 500 624 501 /** 502 * Call the actual INCL routine. 503 */ 504 625 // G4cout <<"Before PNU:" << G4endl; 505 626 // randomGenerator->printSeeds(); 506 G4double jremx, jremy, jremz; 507 pnu(&ibert, &nopart,&izrem,&iarem,&esrem,&erecrem,&alrem,&berem,&garem,&bimpac,&jrem, 508 &jremx, &jremy, &jremz); 627 // Call the actual INCL routine: 628 pnu(&ibert, &nopart,&izrem,&iarem,&esrem,&erecrem,&alrem,&berem,&garem,&bimpac, 629 &jrem, &xjrem, &yjrem, &zjrem); 630 // G4cout <<"After PNU:" << G4endl; 631 // randomGenerator->printSeeds(); 632 G4double mrem = int(zjrem/197.328); // CHECK 633 if (mrem > jrem) mrem=jrem; 634 if (mrem < -jrem) mrem=-jrem; 635 636 // nopart=1; 637 // kind[0]=1; 638 // ep[0]=799.835; 639 // alpha[0]=0.08716; 640 // beta[0]=0.; 641 // gam[0]=0.99619; 642 // izrem=82; 643 // iarem=208; 644 // esrem=200.; 645 // erecrem=0.18870; 646 // alrem=-0.47101; 647 // berem=0.; 648 // garem=0.88213; 649 // bimpac=2.; 509 650 forceAbsor(&nopart, &iarem, &izrem, &esrem, &erecrem, &alrem, &berem, &garem, &jrem); 510 G4double aprf = double(iarem); // mass number of the prefragment511 G4double jprf = 0.0; // angular momentum of the prefragment512 513 // Mean angular momentum of prefragment .514 jprf = 0.165 * std::pow(at,(2.0/3.0)) * aprf * (at - aprf) /(at - 1.0);651 G4double aprf = double(iarem); // mass number of the prefragment 652 G4double jprf = 0.0; // angular momentum of the prefragment 653 654 // Mean angular momentum of prefragment 655 jprf = 0.165 * std::pow(at,(2.0/3.0)) * aprf*(at - aprf)/(at - 1.0); 515 656 if (jprf < 0) { 516 657 jprf = 0.0; 517 658 } 518 659 519 // Reference M. de Jong, Ignatyuk, Schmidt Nuc. Phys A 613, p442, 7th line660 // check m.de jong, ignatyuk, schmidt nuc.phys a 613, pg442, 7th line 520 661 jprf = std::sqrt(2*jprf); 662 521 663 jprf = jrem; 522 varntp->jremn = jrem; // Copy jrem to output tuple. 523 524 G4double numpi = 0; // Number of pions. 525 G4double multn = 0; // Number (multiplicity) of neutrons. 526 G4double multp = 0; // Number (multiplicity) of protons. 527 528 // Ecriture dans le ntuple des particules de cascade (sauf remnant). 529 varntp->ntrack = nopart; // Nombre de particules pour ce tir. 664 varntp->jremn = jrem; // jrem copie dans le ntuple 665 666 G4double numpi = 0; // compteurs de pions, neutrons protons 667 G4double multn = 0; 668 G4double multp = 0; 669 670 // ecriture dans le ntuple des particules de cascade (sauf remnant) 671 varntp->ntrack = nopart; // nombre de particules pour ce tir 672 if(varntp->ntrack >= VARNTPSIZE) { 673 if(verboseLevel > 2) { 674 G4cout <<"G4Incl error: Output data structure not big enough." << G4endl; 675 } 676 } 530 677 varntp->massini = iarem; 531 678 varntp->mzini = izrem; … … 533 680 varntp->bimpact = bimpac; 534 681 535 /** 536 * Three ways to compute the mass of the remnant: 537 * -from the output of the cascade and the canonic mass 538 * -from energy balance (input - all emitted energies) 539 * -following the approximations of the cugnon code (esrem...) 540 */ 541 G4double f0 = calincl->targetA(); 542 G4double f1 = calincl->targetZ(); 543 G4double f2 = calincl->bulletE(); 544 G4double mcorem = mprojo + f2 + abla->pace2(f0, f1) + f0 * uma - f1 * melec; 682 // three ways to compute the mass of the remnant: 683 // -from the output of the cascade and the canonic mass 684 // -from energy balance (input - all emitted energies) 685 // -following the approximations of the cugnon code (esrem...) 686 G4double mcorem = mprojo + calincl->bulletE() + abla->pace2(double(calincl->targetA()),double(calincl->targetZ())) 687 + calincl->targetA()*uma - calincl->targetZ()*melec; 545 688 546 689 G4double pxbil = 0.0; … … 549 692 550 693 if(nopart > -1) { 551 for(G4int j = 0; j < nopart; j++) { 552 varntp->itypcasc[j] = 1; 694 // Fill the projectile spectator variables 695 varntp->masp = ps->a_projspec; 696 varntp->mzsp = ps->z_projspec; 697 varntp->exsp = ps->ex_projspec; 698 varntp->spectatorP1 = ps->p1_projspec; 699 varntp->spectatorP2 = ps->p2_projspec; 700 varntp->spectatorP3 = ps->p3_projspec; 701 varntp->spectatorT = ps->t_projspec; 702 for(G4int j = 0; j <= nopart; j++) { 703 if(ep[j] < 0.0) { // Workaround to avoid negative energies (and taking std::sqrt of a negative number). 704 G4cout <<"G4Incl: Not registering particle with energy: " << ep[j] << G4endl; 705 continue; 706 } 707 if(kind[j] == 0) continue; // Empty particle rows are sometimes produced by lurking indexing problems. We can simply skip those "bad" entries... 708 if(gam[j] > CLHEP::pi) { 709 if(verboseLevel > 2) { 710 G4cout <<"G4Incl: Just avoided floating point exception by using an ugly hack..." << G4endl; 711 } 712 continue; // Avoid floating point exception 713 } 714 715 varntp->itypcasc[j] = 1; // Particle was produced by the cascade 716 // Spectators of composite projectiles (7/2006, AB) 717 // (kind is negative in that case) 718 if(kind[j] <= 0) { // Particle is a projectile spectator that comes directly from the cascade 719 kind[j] *= -1; 720 varntp->itypcasc[j]=-1; 721 // G4cout <<"Spectator registered!" << G4endl; 722 // continue; 723 } 724 553 725 // kind(): 1=proton, 2=neutron, 3=pi+, 4=pi0, 5=pi - 554 726 if(kind[j] == 1) { 555 727 varntp->avv[j] = 1; 556 728 varntp->zvv[j] = 1; 557 varntp->plab[j] = std::sqrt(ep[j] * (ep[j] +1876.5592)); // cugnon729 varntp->plab[j] = std::sqrt(ep[j]*(ep[j]+1876.5592)); // cugnon 558 730 multp = multp + 1; 559 731 mcorem = mcorem - ep[j] - 938.27231; … … 568 740 varntp->zvv[j] = 0; 569 741 varntp->plab[j] = std::sqrt(ep[j]*(ep[j]+1876.5592)); // cugnon 570 //varntp->plab[j] = std::sqrt(ep[j] * (ep[j] + 1879.13126)); // PK mass check571 742 multn = multn + 1; 572 743 mcorem = mcorem - ep[j] - 939.56563; … … 576 747 } 577 748 } 578 749 579 750 if(kind[j] == 3) { 580 751 varntp->avv[j] = -1; … … 656 827 mcorem = mcorem - ep[j] - 3724.818; 657 828 if(verboseLevel > 3) { 658 G4cout <<"G4Incl: He 4produced! " << G4endl;829 G4cout <<"G4Incl: He3 produced! " << G4endl; 659 830 G4cout <<"G4Incl: Momentum: "<< varntp->plab[j] << G4endl; 660 831 } … … 669 840 670 841 if(verboseLevel > 3) { 671 G4cout <<"Momentum: " << varntp->plab[j] << G4endl;672 G4cout <<"Theta: " << varntp->tetlab[j] << G4endl;673 G4cout <<"Phi: " << varntp->philab[j] << G4endl;842 G4cout <<"Momentum: " << varntp->plab[j] << G4endl; 843 G4cout <<"Theta: " << varntp->tetlab[j] << G4endl; 844 G4cout <<"Phi: " << varntp->philab[j] << G4endl; 674 845 } 675 846 } 676 847 677 848 // calcul de la masse (impulsion) du remnant coherente avec la conservation d'energie: 678 pcorem=std::sqrt(erecrem*(erecrem +2.*938.2796*iarem)); // cugnon 679 mcorem = 938.2796*iarem; // cugnon 680 849 pcorem = std::sqrt(erecrem*(erecrem + 2.0 * 938.2796 * iarem)); // cugnon 850 mcorem = 938.2796 * iarem; // cugnon 851 varntp->pcorem = pcorem; 852 varntp->mcorem = mcorem; 681 853 // Note: Il faut negliger l'energie d'excitation (ESREM) pour que le bilan 682 854 // d'impulsion soit correct a la sortie de la cascade.....et prendre la 683 855 // masse MCOREM comme ci-dessus (fausse de ~1GeV par rapport aux tables...) 684 pxrem=pcorem*alrem; 685 pyrem=pcorem*berem; 686 pzrem=pcorem*garem; 687 688 pxbil=pxbil+pxrem; 689 pybil=pybil+pyrem; 690 pzbil=pzbil+pzrem; 691 692 if((std::fabs(pzbil-pbeam) > 5.0) || (std::sqrt(std::pow(pxbil,2)+std::pow(pybil,2)) >= 3.0)) { 856 pxrem = pcorem * alrem; 857 pyrem = pcorem * berem; 858 pzrem = pcorem * garem; 859 varntp->pxrem = pxrem; 860 varntp->pyrem = pyrem; 861 varntp->pzrem = pzrem; 862 pxbil = pxbil + pxrem; 863 pybil = pybil + pyrem; 864 pzbil = pzbil + pzrem; 865 866 // If on purpose, add here the spectator nuclei: 867 if(calincl->bulletType() < 0 && ps->a_projspec != 0) { 868 pxbil=pxbil+ps->p1_projspec; 869 pybil=pybil+ps->p2_projspec; 870 pzbil=pzbil+ps->p3_projspec; 871 } 872 873 if((std::fabs(pzbil - pbeam) > 5.0) || (std::sqrt(std::pow(pxbil,2) + std::pow(pybil,2)) >= 3.0)) { 693 874 if(verboseLevel > 3) { 694 875 G4cout <<"Bad momentum conservation after INCL:" << G4endl; … … 704 885 varntp->iafis = 0; 705 886 706 // varntp->ntrack = varntp->ntrack + 1; // on recopie le remnant dans le ntuple 887 // on recopie le remnant dans le ntuple 888 varntp->ntrack = varntp->ntrack + 1; 707 889 varntp->massini = iarem; 708 890 varntp->mzini = izrem; 709 891 varntp->exini = esrem; 710 varntp->itypcasc[varntp->ntrack] = 1; 711 varntp->avv[varntp->ntrack] = iarem; 712 varntp->zvv[varntp->ntrack]= izrem; 713 varntp->plab[varntp->ntrack] = pcorem; 714 varntp->enerj[varntp->ntrack] = std::sqrt(std::pow(pcorem,2) + std::pow(mcorem,2)) - mcorem; 715 varntp->tetlab[varntp->ntrack] = 180.0*std::acos(garem)/3.141592654; 716 varntp->philab[varntp->ntrack] = 180.0*std::atan2(berem,alrem)/3.141592654; 717 varntp->ntrack++; 718 varntp->mulncasc = varntp->ntrack; 719 varntp->mulnevap = 0; 720 varntp->mulntot = varntp->mulncasc + varntp->mulnevap; 721 if(verboseLevel > 3) { 722 G4cout <<"G4Incl: Returning nucleus fragment. " << G4endl; 723 G4cout <<"G4Incl: Fragment A = " << varntp->avv[varntp->ntrack] << " Z = " << varntp->zvv[varntp->ntrack] << G4endl; 724 G4cout <<"Energy: " << varntp->enerj[varntp->ntrack] << G4endl; 725 G4cout <<"Momentum: " << varntp->plab[varntp->ntrack] << G4endl; 726 G4cout <<"Theta: " << varntp->tetlab[varntp->ntrack] << G4endl; 727 G4cout <<"Phi: " << varntp->philab[varntp->ntrack] << G4endl; 728 } 729 } 730 else { 731 if(nopart == -2) { 732 varntp->ntrack = -2; //FIX: Error flag to remove events containing unphysical events (Ekin > Ebullet). 733 } 734 else { 735 varntp->ntrack = -1; 736 } 892 varntp->pxrem = pxrem; 893 varntp->pyrem = pyrem; 894 varntp->pzrem = pzrem; 895 varntp->mcorem = mcorem; 896 varntp->erecrem = pcorem; 897 varntp->erecrem = erecrem; 898 899 #ifdef G4INCLDEBUG 900 theLogger->fillHistogram1D("bimpact", varntp->bimpact); 901 #endif 902 903 #ifdef G4INCLDEBUG 904 theLogger->fillHistogram1D("mzini", varntp->mzini); 905 #endif 906 } 907 if(nopart == -2) { 908 varntp->ntrack = -2; //FIX: Error flag to remove events containing unphysical events (Ekin > Ebullet). 909 evaporationResult->ntrack = -2; //FIX: Error flag to remove events containing unphysical events (Ekin > Ebullet). 910 } 911 else if(nopart == -1) { 912 varntp->ntrack = -1; 913 evaporationResult->ntrack = -1; 914 } 915 if(verboseLevel > 2) { 916 G4cout << __FILE__ << ":" << __LINE__ << "Dump varntp after combining: " << G4endl; 917 varntp->dump(); 737 918 } 738 919 } … … 8150 8331 } 8151 8332 8152 void G4Incl::standardRandom(G4double *rndm, G4long *seed) 8153 { 8154 (*seed) = (*seed); // Avoid warning during compilation. 8333 void G4Incl::standardRandom(G4double *rndm, G4long*) 8334 { 8155 8335 // Use Geant4 G4UniformRand 8156 8336 // (*rndm) = G4UniformRand(); … … 8550 8730 } // enddo 8551 8731 G4double p_spec2=std::pow(p1_spec,2)+std::pow(p2_spec,2)+std::pow(p3_spec,2); 8552 G4double s_spec = s qrt(std::pow(e_spec,2)-p_spec2);8732 G4double s_spec = std::sqrt(std::pow(e_spec,2)-p_spec2); 8553 8733 8554 8734 // no projectile spectator if a>=4 and a=z or a=n (no dresner breakup) … … 8844 9024 } 8845 9025 9026 // Initialization helper 9027 void G4Incl::clearState() { 9028 G4int epSize = 300; 9029 for(G4int i = 0; i < epSize; ++i) { 9030 kind[i] = 0; 9031 ep[i] = 0.0; 9032 alpha[i] = 0.0; 9033 beta[i] = 0.0; 9034 gam[i] = 0.0; 9035 } 9036 inside_step = 0; 9037 } 9038 8846 9039 // C------------------------------------------------------------------------ 8847 9040 // C Logging … … 8923 9116 // COMMON/BL2/CROIS(19900),K,IND(20000),JND(20000) 8924 9117 // COMMON/BL3/R1,R2,X1(300),X2(300),X3(300),IA1,IA2,RAB2 P-N22270 9118 9119 if(index < 0) { 9120 G4cout <<";; Error! index < 0! index = " << index << G4endl; 9121 return; 9122 } 8925 9123 8926 9124 G4int i_ind1 = bl1->ind1[bl2->ind[index]]; -
trunk/source/processes/hadronic/models/incl/src/G4InclAblaCascadeInterface.cc
r1340 r1347 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclAblaCascadeInterface.cc,v 1. 16 2010/10/29 06:48:43 gunterExp $26 // $Id: G4InclAblaCascadeInterface.cc,v 1.20 2010/11/17 20:19:09 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 116 116 incl->setInput(calincl); 117 117 118 G4InclInput::printProjectileTargetInfo(aTrack, theNucleus);119 calincl->printInfo();118 // G4InclInput::printProjectileTargetInfo(aTrack, theNucleus); 119 // calincl->printInfo(); 120 120 121 121 #ifdef DEBUGINCL … … 439 439 } 440 440 } 441 delete theFermiBreakupResult; 442 theFermiBreakupResult = 0; 443 441 444 if(std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 442 445 G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl; -
trunk/source/processes/hadronic/models/incl/src/G4InclAblaDataFile.cc
r1337 r1347 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclAblaDataFile.cc,v 1. 9 2010/06/14 16:10:01 gcosmoExp $26 // $Id: G4InclAblaDataFile.cc,v 1.10 2010/11/17 20:19:09 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 116 116 vgsldin.close(); 117 117 118 int A = 0, Zbegin = 0, Zend = 0;119 118 G4String str1, str2, str3; 120 119 for(int i = 0; i < 500; i++) { … … 124 123 } 125 124 125 int A = 0, Zbegin = 0, Zend = 0; 126 126 for(int i = 0; i < massnumbers; i++) { 127 127 pace2in >> str1 >> A >> str2 >> Zbegin >> str3 >> Zend; 128 for(int j = Zbegin; j <= Zend; j++) { 129 pace2in >> pace2; 130 setPace2(A, j, pace2); 128 if(Zbegin >= 0 && Zbegin < getPaceCols() && 129 A >= 0 && A < getPaceRows()) { 130 for(int j = Zbegin; j <= Zend; j++) { 131 pace2in >> pace2; 132 setPace2(A, j, pace2); 133 } 131 134 } 132 135 } -
trunk/source/processes/hadronic/models/incl/src/G4InclAblaLightIonInterface.cc
r1340 r1347 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclAblaLightIonInterface.cc,v 1.1 4 2010/10/29 06:48:43 gunterExp $26 // $Id: G4InclAblaLightIonInterface.cc,v 1.16 2010/11/17 20:19:09 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 103 103 } 104 104 105 if(verboseLevel > 1) {106 G4Calincl::printProjectileTargetInfo(aTrack, theNucleus);107 }108 109 105 // Inverse kinematics for targets with Z = 1 and A = 1 110 106 // if(false) { 111 107 G4LorentzRotation toBreit = aTrack.Get4Momentum().boostVector(); 112 108 113 if(theNucleus.GetZ_asInt() == 1 && theNucleus.GetA_asInt() == 1 && G4 Calincl::canUseInverseKinematics(aTrack, theNucleus)) {109 if(theNucleus.GetZ_asInt() == 1 && theNucleus.GetA_asInt() == 1 && G4InclInput::canUseInverseKinematics(aTrack, theNucleus)) { 114 110 G4ParticleDefinition *oldTargetDef = theTableOfParticles->GetIon(theNucleus.GetA_asInt(), theNucleus.GetZ_asInt(), 0.0); 115 111 const G4ParticleDefinition *oldProjectileDef = aTrack.GetDefinition(); 116 112 113 if(oldProjectileDef != 0 && oldTargetDef != 0) { 117 114 G4int oldTargetA = oldTargetDef->GetAtomicMass(); 118 115 G4int newTargetA = oldProjectileDef->GetAtomicMass(); 119 116 G4int newTargetZ = oldProjectileDef->GetAtomicNumber(); 120 117 121 if(newTargetA > 0 && newTargetZ > 0 && oldTargetDef != 0 && oldProjectileDef != 0) {118 if(newTargetA > 0 && newTargetZ > 0) { 122 119 G4Nucleus swappedTarget(oldProjectileDef->GetAtomicMass(), oldProjectileDef->GetAtomicNumber()); 123 120 … … 137 134 G4cout <<"Badly defined target after swapping. Falling back to normal (non-swapped) mode." << G4endl; 138 135 calincl = new G4InclInput(aTrack, theNucleus, false); 136 } 139 137 } 140 138 } else { … … 489 487 } 490 488 } 489 delete theSpectatorFermiBreakupResult; 490 theSpectatorFermiBreakupResult = 0; 491 491 492 if(std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 492 493 G4cout <<"Four-momentum balance after spectator nucleus Fermi break-up:" << G4endl; … … 588 589 } 589 590 } 591 delete theFermiBreakupResult; 592 theFermiBreakupResult = 0; 593 590 594 if(std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 591 595 G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl; … … 697 701 } 698 702 703 delete fermiBreakUp; 699 704 delete calincl; 700 705 calincl = 0; -
trunk/source/processes/hadronic/models/incl/src/G4InclCascadeInterface.cc
r1340 r1347 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclCascadeInterface.cc,v 1.1 2 2010/10/26 02:47:59 kaitanie Exp $26 // $Id: G4InclCascadeInterface.cc,v 1.15 2010/11/17 20:19:09 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 34 34 35 35 #include "G4InclCascadeInterface.hh" 36 #include "G4FermiBreakUp.hh" 36 37 #include "math.h" 37 38 #include "G4GenericIon.hh" 38 39 #include "CLHEP/Random/Random.h" 39 40 40 G4InclCascadeInterface::G4InclCascadeInterface() 41 G4InclCascadeInterface::G4InclCascadeInterface(const G4String& nam) 42 :G4VIntraNuclearTransportModel(nam) 41 43 { 42 44 hazard = new G4Hazard(); … … 45 47 46 48 varntp = new G4VarNtp(); 47 inclInput= 0;49 calincl = 0; 48 50 ws = new G4Ws(); 49 51 mat = new G4Mat(); 50 incl = new G4Incl(hazard, inclInput, ws, mat, varntp); 52 incl = new G4Incl(hazard, calincl, ws, mat, varntp); 53 54 theExcitationHandler = new G4ExcitationHandler; 55 thePrecoModel = new G4PreCompoundModel(theExcitationHandler); 56 57 if(!getenv("G4INCLABLANOFERMIBREAKUP")) { // Use Fermi Break-up by default if it is NOT explicitly disabled 58 incl->setUseFermiBreakUp(true); 59 } 51 60 52 61 verboseLevel = 0; … … 55 64 G4InclCascadeInterface::~G4InclCascadeInterface() 56 65 { 66 delete thePrecoModel; 67 delete theExcitationHandler; 68 57 69 delete hazard; 58 70 delete varntp; 59 delete inclInput;60 71 delete ws; 61 72 delete mat; … … 68 79 69 80 G4int particleI; 81 G4int bulletType = 0; 70 82 71 83 // Print diagnostic messages: 0 = silent, 1 and 2 = verbose … … 91 103 G4DynamicParticle *cascadeParticle = 0; 92 104 G4ParticleDefinition *aParticleDefinition = 0; 105 106 G4ReactionProductVector *thePrecoResult = 0; 107 G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable(); 93 108 94 109 // INCL assumes the projectile particle is going in the direction of … … 102 117 G4LorentzRotation toLabFrame = toZ.inverse(); 103 118 104 inclInput = new G4InclInput(aTrack, theNucleus, false);105 106 119 theResult.Clear(); // Make sure the output data structure is clean. 120 121 calincl = new G4InclInput(aTrack, theNucleus, false); 122 incl->setInput(calincl); 123 124 // G4InclInput::printProjectileTargetInfo(aTrack, theNucleus); 125 // calincl->printInfo(); 107 126 108 127 #ifdef DEBUGINCL … … 117 136 G4double eKinSum = bulletE; 118 137 G4LorentzVector labv = G4LorentzVector(0.0, 0.0, std::sqrt(bulletE*(bulletE + 2.*mass)), bulletE + mass + amass); 138 G4LorentzVector labvA = G4LorentzVector(0.0, 0.0, 0.0, 0.0); 119 139 G4cout <<"Energy in the beginning = " << labv.e() / MeV << G4endl; 120 140 #endif 121 141 122 142 // Check wheter the input is acceptable. 123 if(( inclInput->bulletType() != 0) && ((inclInput->targetA() != 1) && (inclInput->targetZ() != 1))) {143 if((calincl->bulletType() != 0) && ((calincl->targetA() != 1) && (calincl->targetZ() != 1))) { 124 144 ws->nosurf = -2; // Nucleus surface, -2 = Woods-Saxon 125 145 ws->xfoisa = 8; … … 130 150 131 151 mat->nbmat = 1; 132 mat->amat[0] = int( inclInput->targetA());133 mat->zmat[0] = int( inclInput->targetZ());152 mat->amat[0] = int(calincl->targetA()); 153 mat->zmat[0] = int(calincl->targetZ()); 134 154 135 155 incl->initIncl(true); … … 140 160 G4cout <<"G4InclCascadeInterface: Try number = " << nTries << G4endl; 141 161 } 142 incl->processEventIncl( inclInput);162 incl->processEventIncl(calincl); 143 163 144 164 if(verboseLevel > 1) { … … 151 171 * Diagnostic output 152 172 */ 153 G4cout <<"G4InclCascadeInterface: Bullet type: " << inclInput->bulletType() << G4endl;154 G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << inclInput->bulletE() << " MeV" << G4endl;155 156 G4cout <<"G4InclCascadeInterface: Target A: " << inclInput->targetA() << G4endl;157 G4cout <<"G4InclCascadeInterface: Target Z: " << inclInput->targetZ() << G4endl;173 G4cout <<"G4InclCascadeInterface: Bullet type: " << calincl->bulletType() << G4endl; 174 G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl; 175 176 G4cout <<"G4InclCascadeInterface: Target A: " << calincl->targetA() << G4endl; 177 G4cout <<"G4InclCascadeInterface: Target Z: " << calincl->targetZ() << G4endl; 158 178 159 179 if(verboseLevel > 3) { 160 diagdata <<"G4InclCascadeInterface: Bullet type: " << inclInput->bulletType() << G4endl;161 diagdata <<"G4InclCascadeInterface: Bullet energy: " << inclInput->bulletE() << " MeV" << G4endl;180 diagdata <<"G4InclCascadeInterface: Bullet type: " << calincl->bulletType() << G4endl; 181 diagdata <<"G4InclCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl; 162 182 163 diagdata <<"G4InclCascadeInterface: Target A: " << inclInput->targetA() << G4endl; 164 diagdata <<"G4InclCascadeInterface: Target Z: " << inclInput->targetZ() << G4endl; 165 } 166 167 // for(particleI = 0; particleI < varntp->ntrack; particleI++) { 168 // G4cout << n << " " << inclInput->f[6] << " " << inclInput->f[2] << " "; 169 // G4cout << varntp->massini << " " << varntp->mzini << " "; 170 // G4cout << varntp->exini << " " << varntp->mulncasc << " " << varntp->mulnevap << " " << varntp->mulntot << " "; 171 // G4cout << varntp->bimpact << " " << varntp->jremn << " " << varntp->kfis << " " << varntp->estfis << " "; 172 // G4cout << varntp->izfis << " " << varntp->iafis << " " << varntp->ntrack << " " << varntp->itypcasc[particleI] << " "; 173 // G4cout << varntp->avv[particleI] << " " << varntp->zvv[particleI] << " " << varntp->enerj[particleI] << " "; 174 // G4cout << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl; 175 // // For diagnostic output 176 // if(verboseLevel > 3) { 177 // diagdata << n << " " << inclInput->f[6] << " " << inclInput->f[2] << " "; 178 // diagdata << varntp->massini << " " << varntp->mzini << " "; 179 // diagdata << varntp->exini << " " << varntp->mulncasc << " " << varntp->mulnevap << " " << varntp->mulntot << " "; 180 // diagdata << varntp->bimpact << " " << varntp->jremn << " " << varntp->kfis << " " << varntp->estfis << " "; 181 // diagdata << varntp->izfis << " " << varntp->iafis << " " << varntp->ntrack << " "; 182 // diagdata << varntp->itypcasc[particleI] << " "; 183 // diagdata << varntp->avv[particleI] << " " << varntp->zvv[particleI] << " " << varntp->enerj[particleI] << " "; 184 // diagdata << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl; 185 // } 186 // } 183 diagdata <<"G4InclCascadeInterface: Target A: " << calincl->targetA() << G4endl; 184 diagdata <<"G4InclCascadeInterface: Target Z: " << calincl->targetZ() << G4endl; 185 } 186 187 187 } 188 188 … … 196 196 197 197 theResult.SetStatusChange(stopAndKill); 198 aParticleDefinition = G4InclInput::getParticleDefinition(inclInput->bulletType()); 199 200 cascadeParticle = new G4DynamicParticle(); 201 cascadeParticle->SetDefinition(aParticleDefinition); 202 cascadeParticle->Set4Momentum(aTrack.Get4Momentum()); 203 theResult.AddSecondary(cascadeParticle); 198 199 G4int bulletType = calincl->bulletType(); 200 aParticleDefinition = G4InclInput::getParticleDefinition(bulletType); 201 202 if(aParticleDefinition != 0) { 203 cascadeParticle = new G4DynamicParticle(); 204 cascadeParticle->SetDefinition(aParticleDefinition); 205 cascadeParticle->Set4Momentum(aTrack.Get4Momentum()); 206 theResult.AddSecondary(cascadeParticle); 207 } 204 208 } 205 209 … … 209 213 210 214 #ifdef DEBUGINCL 211 G4cout << "E [MeV]" << std::setw(12) << " Ekin [MeV]" << std::setw(12) << " E* [MeV]" << std::setw(12) << "Px [MeV]" << std::setw(12) << " Py [MeV]" << std::setw(12) << "Pz [MeV]" << std::setw(12) << "Pt [MeV]" << std::setw(12) << "A" << std::setw(12) << "Z" << G4endl; 215 G4cout << "E [MeV]" << std::setw(12) 216 << " Ekin [MeV]" << std::setw(12) 217 << "Px [MeV]" << std::setw(12) 218 << " Py [MeV]" << std::setw(12) 219 << "Pz [MeV]" << std::setw(12) 220 << "Pt [MeV]" << std::setw(12) 221 << "A" << std::setw(12) 222 << "Z" << G4endl; 212 223 #endif 213 224 … … 276 287 if((varntp->avv[particleI] > 1) && (varntp->zvv[particleI] >= 1)) { // Nucleus fragment 277 288 G4ParticleDefinition * aIonDef = 0; 278 G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();279 289 280 290 G4int A = G4int(varntp->avv[particleI]); … … 315 325 G4double p = mom.mag(); 316 326 labv -= fm; 317 G4double px = mom.x() * MeV; 318 G4double py = mom.y() * MeV; 319 G4double pz = mom.z() * MeV; 327 if(varntp->avv[particleI] > 1) { 328 labvA += fm; 329 } 330 G4double px = mom.x() * MeV; 331 G4double py = mom.y() * MeV; 332 G4double pz = mom.z() * MeV; 320 333 G4double pt = std::sqrt(px*px+py*py); 321 334 G4double e = fm.e(); … … 330 343 G4cout << fm.e() / MeV 331 344 << std::setw(12) << cascadeParticle->GetKineticEnergy() / MeV 332 << std::setw(12) << exE / MeV333 345 << std::setw(12) << mom.x() / MeV 334 346 << std::setw(12) << mom.y() / MeV … … 351 363 } 352 364 } 365 366 G4double nuclearMass = G4NucleiProperties::GetNuclearMass(G4int(varntp->massini), G4int(varntp->mzini)) + varntp->exini * MeV; 367 G4LorentzVector fragmentMomentum(varntp->pxrem * MeV, varntp->pyrem * MeV, varntp->pzrem * MeV, 368 varntp->erecrem * MeV + nuclearMass); 369 G4double momentumScaling = G4InclUtils::calculate4MomentumScaling(G4int(varntp->massini), G4int(varntp->mzini), 370 varntp->exini, 371 varntp->erecrem, 372 varntp->pxrem, 373 varntp->pyrem, 374 varntp->pzrem); 375 G4LorentzVector p4(momentumScaling * varntp->pxrem * MeV, momentumScaling * varntp->pyrem * MeV, 376 momentumScaling * varntp->pzrem * MeV, 377 varntp->erecrem + nuclearMass); 378 379 // For four-momentum, baryon number and charge conservation check: 380 G4LorentzVector fourMomentumBalance = p4; 381 G4int baryonNumberBalance = G4int(varntp->massini); 382 G4int chargeBalance = G4int(varntp->mzini); 383 384 G4LorentzRotation toFragmentZ; 385 toFragmentZ.rotateZ(-p4.theta()); 386 toFragmentZ.rotateY(-p4.phi()); 387 G4LorentzRotation toFragmentLab = toFragmentZ.inverse(); 388 p4 *= toFragmentZ; 389 390 G4LorentzVector p4rest = p4; 391 p4rest.boost(-p4.boostVector()); 392 if(verboseLevel > 0) { 393 G4cout <<"Cascade remnant nucleus:" << G4endl; 394 G4cout <<"p4: " << G4endl; 395 G4cout <<" px: " << p4.px() <<" py: " << p4.py() <<" pz: " << p4.pz() << G4endl; 396 G4cout <<" E = " << p4.e() << G4endl; 397 398 G4cout <<"p4rest: " << G4endl; 399 G4cout <<" px: " << p4rest.px() <<" py: " << p4rest.py() <<" pz: " << p4rest.pz() << G4endl; 400 G4cout <<" E = " << p4rest.e() << G4endl; 401 } 402 403 G4Fragment theCascadeRemnant(G4int(varntp->massini), G4int(varntp->mzini), p4rest); 404 thePrecoResult = thePrecoModel->DeExcite(theCascadeRemnant); 405 if(thePrecoResult != 0) { 406 G4ReactionProductVector::iterator fragment; 407 for(fragment = thePrecoResult->begin(); fragment != thePrecoResult->end(); fragment++) { 408 G4ParticleDefinition *theFragmentDefinition = (*fragment)->GetDefinition(); 409 410 if(theFragmentDefinition != 0) { 411 G4DynamicParticle *theFragment = new G4DynamicParticle(theFragmentDefinition, (*fragment)->GetMomentum()); 412 G4LorentzVector labMomentum = theFragment->Get4Momentum(); 413 labMomentum.boost(p4.boostVector()); 414 labMomentum *= toFragmentLab; 415 labMomentum *= toLabFrame; 416 theFragment->Set4Momentum(labMomentum); 417 fourMomentumBalance -= theFragment->Get4Momentum(); 418 baryonNumberBalance -= theFragmentDefinition->GetAtomicMass(); 419 chargeBalance -= theFragmentDefinition->GetAtomicNumber(); 420 if(verboseLevel > 0) { 421 G4cout <<"Resulting fragment: " << G4endl; 422 G4cout <<" kinetic energy = " << theFragment->GetKineticEnergy() / MeV << " MeV" << G4endl; 423 G4cout <<" momentum = " << theFragment->GetMomentum().mag() / MeV << " MeV" << G4endl; 424 } 425 theResult.AddSecondary(theFragment); 426 } else { 427 G4cout <<"G4InclCascadeInterface: Error. Fragment produced by Fermi break-up does not exist." << G4endl; 428 G4cout <<"Resulting fragment: " << G4endl; 429 G4cout <<" momentum = " << (*fragment)->GetMomentum().mag() / MeV << " MeV" << G4endl; 430 } 431 } 432 delete thePrecoResult; 433 thePrecoResult = 0; 434 435 if(verboseLevel > 1 && std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 436 G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl; 437 G4cout <<"Magnitude: " << fourMomentumBalance.mag() / MeV << " MeV" << G4endl; 438 G4cout <<"Vector components (px, py, pz, E) = (" 439 << fourMomentumBalance.px() << ", " 440 << fourMomentumBalance.py() << ", " 441 << fourMomentumBalance.pz() << ", " 442 << fourMomentumBalance.e() << ")" << G4endl; 443 } 444 if(baryonNumberBalance != 0 && verboseLevel > 1) { 445 G4cout <<"Baryon number balance after remnant nucleus Fermi break-up: " << baryonNumberBalance << G4endl; 446 } 447 if(chargeBalance != 0 && verboseLevel > 1) { 448 G4cout <<"Charge balance after remnant nucleus Fermi break-up: " << chargeBalance << G4endl; 449 } 450 } 451 // } // if(needsFermiBreakUp) 452 353 453 #ifdef DEBUGINCL 354 454 G4cout <<"--------------------------------------------------------------------------------" << G4endl; 355 455 G4double pt = std::sqrt(std::pow(labv.x(), 2) + std::pow(labv.y(), 2)); 356 G4cout << labv.e() / MeV << std::setw(12) << eKinSum / MeV << std::setw(12) << labv.x() << std::setw(12) << labv.y() << std::setw(12) << labv.z() << std::setw(12) << pt / MeV << std::setw(12) << baryonNumber << std::setw(12) << chargeNumber << " totals" << G4endl; 456 G4double ptA = std::sqrt(std::pow(labvA.x(), 2) + std::pow(labvA.y(), 2)); 457 G4cout << labv.e() / MeV << std::setw(12) 458 << eKinSum / MeV << std::setw(12) 459 << labv.x() / MeV << std::setw(12) 460 << labv.y() / MeV << std::setw(12) 461 << labv.z() / MeV << std::setw(12) 462 << pt / MeV << std::setw(12) 463 << baryonNumber << std::setw(12) 464 << chargeNumber << " totals" << G4endl; 465 G4cout << " - " << std::setw(12) 466 << " - " << std::setw(12) 467 << labvA.x() / MeV << std::setw(12) 468 << labvA.y() / MeV << std::setw(12) 469 << labvA.z() / MeV << std::setw(12) 470 << ptA / MeV << std::setw(12) 471 << " - " << std::setw(12) << " - " << " totals ABLA" << G4endl; 357 472 G4cout << G4endl; 358 473 359 474 if(verboseLevel > 3) { 360 475 if(baryonNumber != 0) { … … 369 484 } 370 485 #endif 371 486 372 487 varntp->ntrack = 0; // Clean up the number of generated particles in the event. 373 488 } … … 379 494 theResult.SetStatusChange(stopAndKill); 380 495 381 G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();382 496 cascadeParticle = new G4DynamicParticle(theTableOfParticles->FindParticle(aTrack.GetDefinition()), aTrack.Get4Momentum()); 383 497 … … 388 502 } 389 503 if(verboseLevel > 3) { 390 diagdata <<"ERROR G4InclCascadeInterface: Processing event number (internal) failed " << eventNumber << G4endl;391 } 392 393 if( inclInput->bulletType()== 0) {504 diagdata <<"ERROR G4InclCascadeInterface: Error processing event number (internal) failed " << eventNumber << G4endl; 505 } 506 507 if(bulletType == 0) { 394 508 if(verboseLevel > 1) { 395 509 G4cout <<"G4InclCascadeInterface: Unknown bullet type" << G4endl; … … 402 516 } 403 517 404 if(( inclInput->targetA() == 1) && (inclInput->targetZ() == 1)) { // Unsupported target518 if((calincl->targetA() == 1) && (calincl->targetZ() == 1)) { // Unsupported target 405 519 if(verboseLevel > 1) { 406 520 G4cout <<"Unsupported target: " << G4endl; 407 G4cout <<"Target A: " << inclInput->targetA() << G4endl;408 G4cout <<"TargetZ: " << inclInput->targetZ() << G4endl;521 G4cout <<"Target A: " << calincl->targetA() << G4endl; 522 G4cout <<"TargetZ: " << calincl->targetZ() << G4endl; 409 523 } 410 524 if(verboseLevel > 3) { 411 525 diagdata <<"Unsupported target: " << G4endl; 412 diagdata <<"Target A: " << inclInput->targetA() << G4endl;413 diagdata <<"TargetZ: " << inclInput->targetZ() << G4endl;414 } 415 } 416 417 if( inclInput->bulletE() < 100) { // INCL does not support E < 100 MeV.526 diagdata <<"Target A: " << calincl->targetA() << G4endl; 527 diagdata <<"TargetZ: " << calincl->targetZ() << G4endl; 528 } 529 } 530 531 if(calincl->bulletE() < 100) { // INCL does not support E < 100 MeV. 418 532 if(verboseLevel > 1) { 419 G4cout <<"Unsupported bullet energy: " << inclInput->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;533 G4cout <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl; 420 534 G4cout <<"WARNING: Returning the original bullet with original energy back to Geant4." << G4endl; 421 535 } 422 536 if(verboseLevel > 3) { 423 diagdata <<"Unsupported bullet energy: " << inclInput->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;537 diagdata <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl; 424 538 } 425 539 } … … 430 544 } 431 545 546 delete calincl; 547 calincl = 0; 432 548 return &theResult; 433 549 } -
trunk/source/processes/hadronic/models/incl/src/G4InclLightIonInterface.cc
r1340 r1347 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclLightIonInterface.cc,v 1.1 2 2010/10/26 02:47:59 kaitanie Exp $26 // $Id: G4InclLightIonInterface.cc,v 1.15 2010/11/17 20:19:09 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 31 31 // Aatos Heikkinen, HIP (project coordination) 32 32 33 #include <vector> 34 33 35 #include "G4InclLightIonInterface.hh" 36 #include "G4FermiBreakUp.hh" 34 37 #include "math.h" 35 38 #include "G4GenericIon.hh" … … 39 42 { 40 43 hazard = new G4Hazard(); 44 41 45 const G4long* table_entry = CLHEP::HepRandom::getTheSeeds(); // Get random seed from CLHEP. 42 46 hazard->ial = (*table_entry); 47 48 theExcitationHandler = new G4ExcitationHandler; 49 thePrecoModel = new G4PreCompoundModel(theExcitationHandler); 43 50 44 51 varntp = new G4VarNtp(); … … 47 54 mat = new G4Mat(); 48 55 incl = new G4Incl(hazard, calincl, ws, mat, varntp); 49 56 useProjectileSpectator = true; 57 useFermiBreakup = true; 58 incl->setUseProjectileSpectators(useProjectileSpectator); 59 if(!getenv("G4INCLABLANOFERMIBREAKUP")) { // Use Fermi Break-up by default if it is NOT explicitly disabled 60 incl->setUseFermiBreakUp(true); 61 useFermiBreakup = true; 62 } 50 63 verboseLevel = 0; 64 if(getenv("G4INCLVERBOSE")) { 65 verboseLevel = 1; 66 } 51 67 } 52 68 53 69 G4InclLightIonInterface::~G4InclLightIonInterface() 54 70 { 71 delete thePrecoModel; 72 delete theExcitationHandler; 73 55 74 delete hazard; 56 75 delete varntp; … … 63 82 G4HadFinalState* G4InclLightIonInterface::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& theNucleus) 64 83 { 84 // const G4bool useFermiBreakup = false; 65 85 G4int maxTries = 200; 66 86 67 G4int bulletType = 0;68 87 G4int particleI; 69 88 70 // Print diagnostic messages: 0 = silent, 1 and 2 = verbose 71 verboseLevel = 0; 89 G4int baryonNumberBalanceInINCL = 0; 90 G4int chargeNumberBalanceInINCL = 0; 91 92 G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable(); 72 93 73 94 // Increase the event number: 74 95 eventNumber++; 75 96 97 // Clean up the INCL input 98 if(calincl != 0) { 99 delete calincl; 100 calincl = 0; 101 } 102 76 103 if (verboseLevel > 1) { 77 104 G4cout << " >>> G4InclLightIonInterface::ApplyYourself called" << G4endl; … … 82 109 } 83 110 84 // INCL4 needs the energy in units MeV 85 G4double bulletE = aTrack.GetKineticEnergy() / MeV; 86 87 G4int targetA = theNucleus.GetA_asInt(); 88 G4int targetZ = theNucleus.GetZ_asInt(); 111 // Inverse kinematics for targets with Z = 1 and A = 1 112 // if(false) { 113 G4LorentzRotation toBreit = aTrack.Get4Momentum().boostVector(); 114 115 if(theNucleus.GetZ_asInt() == 1 && theNucleus.GetA_asInt() == 1 && G4InclInput::canUseInverseKinematics(aTrack, theNucleus)) { 116 G4ParticleDefinition *oldTargetDef = theTableOfParticles->GetIon(theNucleus.GetA_asInt(), theNucleus.GetZ_asInt(), 0.0); 117 const G4ParticleDefinition *oldProjectileDef = aTrack.GetDefinition(); 118 119 if(oldTargetDef != 0 && oldProjectileDef != 0) { 120 G4int oldTargetA = oldTargetDef->GetAtomicMass(); 121 G4int newTargetA = oldProjectileDef->GetAtomicMass(); 122 G4int newTargetZ = oldProjectileDef->GetAtomicNumber(); 123 124 if(newTargetA > 0 && newTargetZ > 0) { 125 G4Nucleus swappedTarget(oldProjectileDef->GetAtomicMass(), oldProjectileDef->GetAtomicNumber()); 126 127 // G4cout <<"Original projectile kinE = " << aTrack.GetKineticEnergy() / MeV << G4endl; 128 129 // We need the same energy/nucleon. 130 G4double projectileE = ((aTrack.GetKineticEnergy() / MeV) / newTargetA) * oldTargetA * MeV; 131 132 // G4cout <<"projectileE = " << projectileE << G4endl; 133 G4DynamicParticle swappedProjectileParticle(oldTargetDef, G4ThreeVector(0.0, 0.0, 1.0), projectileE); 134 const G4LorentzVector swapped4Momentum = (swappedProjectileParticle.Get4Momentum()*=toBreit); 135 swappedProjectileParticle.Set4Momentum(swapped4Momentum); 136 const G4HadProjectile swappedProjectile(swappedProjectileParticle); 137 // G4cout <<"New projectile kinE = " << swappedProjectile.GetKineticEnergy() / MeV << G4endl; 138 calincl = new G4InclInput(swappedProjectile, swappedTarget, true); 139 } else { 140 G4cout <<"Badly defined target after swapping. Falling back to normal (non-swapped) mode." << G4endl; 141 calincl = new G4InclInput(aTrack, theNucleus, false); 142 } 143 } 144 } else { 145 calincl = new G4InclInput(aTrack, theNucleus, false); 146 } 89 147 90 148 G4double eKin; … … 103 161 G4LorentzRotation toLabFrame = toZ.inverse(); 104 162 163 /* 164 G4cout <<"Projectile theta = " << projectileMomentum.theta() << " phi = " << projectileMomentum.phi() << G4endl; 165 G4cout <<"Projectile momentum " 166 << "(px = " << projectileMomentum.px() 167 << ", py = " << projectileMomentum.py() 168 << ", pz = " << projectileMomentum.pz() << ")" << G4endl; 169 G4cout << "Projectile energy = " << bulletE << " MeV" << G4endl; 170 */ 171 172 G4ReactionProductVector *thePrecoResult = 0; 173 G4ReactionProductVector *theSpectatorPrecoResult = 0; 174 105 175 theResult.Clear(); // Make sure the output data structure is clean. 176 177 std::vector<G4DynamicParticle*> result; // Temporary list for the results 106 178 107 179 // Map Geant4 particle types to corresponding INCL4 types. 108 180 enum bulletParticleType {nucleus = 0, proton = 1, neutron = 2, pionPlus = 3, pionZero = 4, 109 pionMinus = 5, deuteron = 6, triton = 7, he3 = 8, he4 = 9}; 110 111 // Coding particles for use with INCL4 and ABLA 112 if (aTrack.GetDefinition() == G4Deuteron::Deuteron() ) bulletType = deuteron; 113 if (aTrack.GetDefinition() == G4Triton::Triton() ) bulletType = triton; 114 if (aTrack.GetDefinition() == G4He3::He3() ) bulletType = he3; 115 if (aTrack.GetDefinition() == G4Alpha::Alpha() ) bulletType = he4; 116 117 calincl = new G4InclInput(); 181 pionMinus = 5, deuteron = 6, triton = 7, he3 = 8, he4 = 9, 182 c12 = -12}; // Carbon beam support. 183 184 G4int bulletType = calincl->bulletType(); 185 chargeNumberBalanceInINCL = calincl->targetZ(); 186 baryonNumberBalanceInINCL = calincl->targetA(); 187 188 // G4cout <<"Type of the projectile (INCL projectile code): " << bulletType << G4endl; 189 190 if(bulletType == proton) { 191 chargeNumberBalanceInINCL += 1; 192 baryonNumberBalanceInINCL += 1; 193 } else if(bulletType == neutron) { 194 baryonNumberBalanceInINCL += 1; 195 } else if(bulletType == pionPlus) { //Note: positive pion doesn't contribute to the baryon and charge number counters 196 chargeNumberBalanceInINCL += 1; 197 } else if(bulletType == pionMinus) { 198 chargeNumberBalanceInINCL -= 1; 199 } else if(bulletType == deuteron) { 200 chargeNumberBalanceInINCL += 1; 201 baryonNumberBalanceInINCL += 2; 202 } else if(bulletType == triton) { 203 chargeNumberBalanceInINCL += 1; 204 baryonNumberBalanceInINCL += 3; 205 } else if(bulletType == he3) { 206 chargeNumberBalanceInINCL += 2; 207 baryonNumberBalanceInINCL += 3; 208 } else if(bulletType == he4) { 209 chargeNumberBalanceInINCL += 2; 210 baryonNumberBalanceInINCL += 4; 211 } if(bulletType == c12) { 212 chargeNumberBalanceInINCL += 6; 213 baryonNumberBalanceInINCL += 12; 214 } if(bulletType == -666) { 215 chargeNumberBalanceInINCL += calincl->extendedProjectileZ(); 216 baryonNumberBalanceInINCL += calincl->extendedProjectileA(); 217 } 118 218 119 219 // Check wheter the input is acceptable. 120 if((bulletType != 0) && (( targetA != 1) && (targetZ!= 1))) {220 if((bulletType != 0) && ((calincl->targetA() != 1) && (calincl->targetZ() != 1))) { 121 221 ws->nosurf = -2; // Nucleus surface, -2 = Woods-Saxon 122 222 ws->xfoisa = 8; … … 128 228 mat->nbmat = 1; 129 229 mat->amat[0] = int(calincl->targetA()); 130 mat->zmat[0] = int(calincl->targetZ()); 131 230 mat->zmat[0] = int(calincl->targetA()); 231 232 incl->setInput(calincl); 132 233 incl->initIncl(true); 133 234 … … 148 249 * Diagnostic output 149 250 */ 150 G4cout <<"G4InclLightIonInterface: Bullet type: " << bulletType << G4endl; 151 G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << bulletE << " MeV" << G4endl; 152 153 G4cout <<"G4InclLightIonInterface: Target A: " << targetA << G4endl; 154 G4cout <<"G4InclLightIonInterface: Target Z: " << targetZ << G4endl; 251 G4cout <<"G4InclLightIonInterface: Bullet type: " << calincl->bulletType() << G4endl; 252 G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl; 253 if(bulletType == -666) { 254 G4cout <<" Extended projectile: A = " << calincl->extendedProjectileA() 255 <<" Z = " << calincl->extendedProjectileZ() << G4endl; 256 } 257 258 G4cout <<"G4InclLightIonInterface: Target A: " << calincl->targetA() << G4endl; 259 G4cout <<"G4InclLightIonInterface: Target Z: " << calincl->targetZ() << G4endl; 155 260 156 261 if(verboseLevel > 3) { 157 diagdata <<"G4InclLightIonInterface: Bullet type: " << bulletType<< G4endl;158 diagdata <<"G4InclLightIonInterface: Bullet energy: " << bulletE<< " MeV" << G4endl;262 diagdata <<"G4InclLightIonInterface: Bullet type: " << calincl->bulletType() << G4endl; 263 diagdata <<"G4InclLightIonInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl; 159 264 160 diagdata <<"G4InclLightIonInterface: Target A: " << targetA<< G4endl;161 diagdata <<"G4InclLightIonInterface: Target Z: " << targetZ<< G4endl;265 diagdata <<"G4InclLightIonInterface: Target A: " << calincl->targetA() << G4endl; 266 diagdata <<"G4InclLightIonInterface: Target Z: " << calincl->targetZ() << G4endl; 162 267 } 163 268 } … … 175 280 if(bulletType == proton) { 176 281 aParticleDefinition = G4Proton::ProtonDefinition(); 177 } 178 if(bulletType == neutron) { 282 } else if(bulletType == neutron) { 179 283 aParticleDefinition = G4Neutron::NeutronDefinition(); 180 } 181 if(bulletType == pionPlus) { 284 } else if(bulletType == pionPlus) { 182 285 aParticleDefinition = G4PionPlus::PionPlusDefinition(); 183 } 184 if(bulletType == pionZero) { 286 } else if(bulletType == pionZero) { 185 287 aParticleDefinition = G4PionZero::PionZeroDefinition(); 186 } 187 if(bulletType == pionMinus) { 288 } else if(bulletType == pionMinus) { 188 289 aParticleDefinition = G4PionMinus::PionMinusDefinition(); 189 } 190 191 cascadeParticle = new G4DynamicParticle(); 192 cascadeParticle->SetDefinition(aParticleDefinition); 193 cascadeParticle->Set4Momentum(aTrack.Get4Momentum()); 194 theResult.AddSecondary(cascadeParticle); 290 } else if(bulletType == deuteron) { 291 aParticleDefinition = G4Deuteron::DeuteronDefinition(); 292 } else if(bulletType == triton) { 293 aParticleDefinition = G4Triton::TritonDefinition(); 294 } else if(bulletType == he3) { 295 aParticleDefinition = G4He3::He3Definition(); 296 } else if(bulletType == he4) { 297 aParticleDefinition = G4Alpha::AlphaDefinition(); 298 } else { // Particle was not recognized. Probably an unsupported particle was given as input 299 aParticleDefinition = 0; 300 } 301 302 if(aParticleDefinition != 0) { 303 cascadeParticle = new G4DynamicParticle(); 304 cascadeParticle->SetDefinition(aParticleDefinition); 305 cascadeParticle->Set4Momentum(aTrack.Get4Momentum()); 306 result.push_back(cascadeParticle); 307 } 195 308 } 196 309 … … 199 312 theResult.SetStatusChange(stopAndKill); 200 313 201 for(particleI = 0; particleI < varntp->ntrack; particleI++) { // Loop through the INCL4+ABLA output.314 for(particleI = 0; particleI <= varntp->ntrack; particleI++) { // Loop through the INCL4+ABLA output. 202 315 // Get energy/momentum and construct momentum vector in INCL4 coordinates. 316 // if(varntp->itypcasc[particleI] == -1) continue; // Avoid nucleons that are part of the spectator 317 if(varntp->avv[particleI] == 0 && varntp->zvv[particleI] == 0) continue; 203 318 momx = varntp->plab[particleI]*std::sin(varntp->tetlab[particleI]*CLHEP::pi/180.0)*std::cos(varntp->philab[particleI]*CLHEP::pi/180.0)*MeV; 204 319 momy = varntp->plab[particleI]*std::sin(varntp->tetlab[particleI]*CLHEP::pi/180.0)*std::sin(varntp->philab[particleI]*CLHEP::pi/180.0)*MeV; … … 222 337 new G4DynamicParticle(G4Proton::ProtonDefinition(), momDirection, eKin); 223 338 particleIdentified++; 339 baryonNumberBalanceInINCL -= 1; 340 chargeNumberBalanceInINCL -= 1; 224 341 } 225 342 … … 228 345 new G4DynamicParticle(G4Neutron::NeutronDefinition(), momDirection, eKin); 229 346 particleIdentified++; 347 baryonNumberBalanceInINCL -= 1; 230 348 } 231 349 … … 234 352 new G4DynamicParticle(G4PionPlus::PionPlusDefinition(), momDirection, eKin); 235 353 particleIdentified++; 354 chargeNumberBalanceInINCL -= 1; 236 355 } 237 356 … … 240 359 new G4DynamicParticle(G4PionZero::PionZeroDefinition(), momDirection, eKin); 241 360 particleIdentified++; 361 chargeNumberBalanceInINCL -= 0; 242 362 } 243 363 … … 246 366 new G4DynamicParticle(G4PionMinus::PionMinusDefinition(), momDirection, eKin); 247 367 particleIdentified++; 368 chargeNumberBalanceInINCL -= -1; 248 369 } 249 370 250 371 if((varntp->avv[particleI] > 1) && (varntp->zvv[particleI] >= 1)) { // Nucleus fragment 251 G4ParticleDefinition * aIonDef = 0; 252 G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable(); 372 G4ParticleDefinition * aIonDef = 0; 253 373 254 374 G4int A = G4int(varntp->avv[particleI]); … … 273 393 new G4DynamicParticle(aIonDef, momDirection, eKin); 274 394 particleIdentified++; 395 baryonNumberBalanceInINCL -= A; 396 chargeNumberBalanceInINCL -= Z; 275 397 } 276 398 } 277 399 278 400 if(particleIdentified == 1) { // Particle identified properly. 279 280 theResult.AddSecondary(cascadeParticle); // Put data into G4HadFinalState.401 cascadeParticle->Set4Momentum(cascadeParticle->Get4Momentum()*=toLabFrame); 402 result.push_back(cascadeParticle); 281 403 } 282 404 else { // Particle identification failed. … … 292 414 } 293 415 416 // Spectator nucleus Fermi break-up 417 if(useFermiBreakup && useProjectileSpectator && varntp->masp > 1) { 418 baryonNumberBalanceInINCL -= G4int(varntp->masp); 419 G4double nuclearMass = G4NucleiProperties::GetNuclearMass(G4int(varntp->masp), G4int(varntp->mzsp)) + varntp->exsp * MeV; 420 // Use momentum scaling to compensate for different masses in G4 and INCL: 421 G4double momentumScaling = G4InclUtils::calculate4MomentumScaling(G4int(varntp->masp), 422 G4int(varntp->mzsp), 423 varntp->exsp, 424 varntp->spectatorT, 425 varntp->spectatorP1, 426 varntp->spectatorP2, 427 varntp->spectatorP3); 428 G4LorentzVector p4(momentumScaling * varntp->spectatorP1 * MeV, momentumScaling * varntp->spectatorP2 * MeV, 429 momentumScaling * varntp->spectatorP3 * MeV, 430 varntp->spectatorT * MeV + nuclearMass); 431 // Four-momentum, baryon number and charge balance: 432 G4LorentzVector fourMomentumBalance = p4; 433 G4int baryonNumberBalance = G4int(varntp->masp); 434 chargeNumberBalanceInINCL -= G4int(varntp->mzsp); 435 G4int chargeBalance = G4int(varntp->mzsp); 436 437 G4LorentzRotation toFragmentZ; 438 // Assume that Fermi breakup uses Z as the direction of the projectile 439 toFragmentZ.rotateZ(-p4.theta()); 440 toFragmentZ.rotateY(-p4.phi()); 441 G4LorentzRotation toFragmentLab = toFragmentZ.inverse(); 442 p4 *= toFragmentZ; 443 444 G4LorentzVector p4rest = p4; 445 p4rest.boost(-p4.boostVector()); 446 if(verboseLevel > 0) { 447 G4cout <<"Spectator nucleus:" << G4endl; 448 G4cout <<"p4: " << G4endl; 449 G4cout <<" px: " << p4.px() <<" py: " << p4.py() <<" pz: " << p4.pz() << G4endl; 450 G4cout <<" E = " << p4.e() << G4endl; 451 G4cout <<"p4rest: " << G4endl; 452 G4cout <<" px: " << p4rest.px() <<" py: " << p4rest.py() <<" pz: " << p4rest.pz() << G4endl; 453 G4cout <<" E = " << p4rest.e() << G4endl; 454 } 455 G4Fragment theSpectatorNucleus(G4int(varntp->masp), G4int(varntp->mzsp), p4rest); 456 theSpectatorPrecoResult = thePrecoModel->DeExcite(theSpectatorNucleus); 457 if(theSpectatorPrecoResult != 0) { 458 G4ReactionProductVector::iterator fragment; 459 for(fragment = theSpectatorPrecoResult->begin(); fragment != theSpectatorPrecoResult->end(); fragment++) { 460 G4ParticleDefinition *theFragmentDefinition = (*fragment)->GetDefinition(); 461 462 if(theFragmentDefinition != 0) { 463 G4DynamicParticle *theFragment = new G4DynamicParticle(theFragmentDefinition, (*fragment)->GetMomentum()); 464 G4LorentzVector labMomentum = theFragment->Get4Momentum(); 465 labMomentum.boost(p4.boostVector()); 466 labMomentum *= toFragmentLab; 467 labMomentum *= toLabFrame; 468 theFragment->Set4Momentum(labMomentum); 469 fourMomentumBalance -= theFragment->Get4Momentum(); 470 baryonNumberBalance -= theFragmentDefinition->GetAtomicMass(); 471 chargeBalance -= theFragmentDefinition->GetAtomicNumber(); 472 if(verboseLevel > 0) { 473 G4cout <<"Resulting fragment: " << G4endl; 474 G4cout <<" kinetic energy = " << theFragment->GetKineticEnergy() / MeV << " MeV" << G4endl; 475 G4cout <<" momentum = " << theFragment->GetMomentum().mag() / MeV << " MeV" << G4endl; 476 } 477 theResult.AddSecondary(theFragment); 478 } else { 479 G4cout <<"G4InclCascadeInterface: Error. Fragment produced by Fermi break-up does not exist." << G4endl; 480 G4cout <<"Resulting fragment: " << G4endl; 481 G4cout <<" momentum = " << (*fragment)->GetMomentum().mag() / MeV << " MeV" << G4endl; 482 } 483 } 484 delete theSpectatorPrecoResult; 485 theSpectatorPrecoResult = 0; 486 487 if(verboseLevel > 1 && std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 488 G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl; 489 G4cout <<"Magnitude: " << fourMomentumBalance.mag() / MeV << " MeV" << G4endl; 490 G4cout <<"Vector components (px, py, pz, E) = (" 491 << fourMomentumBalance.px() << ", " 492 << fourMomentumBalance.py() << ", " 493 << fourMomentumBalance.pz() << ", " 494 << fourMomentumBalance.e() << ")" << G4endl; 495 } 496 if(baryonNumberBalance != 0 && verboseLevel > 1) { 497 G4cout <<"Baryon number balance after remnant nucleus Fermi break-up: " << baryonNumberBalance << G4endl; 498 } 499 if(chargeBalance != 0 && verboseLevel > 1) { 500 G4cout <<"Charge balance after remnant nucleus Fermi break-up: " << chargeBalance << G4endl; 501 } 502 } 503 } 504 505 // Finally do Fermi break-up if needed 506 if(varntp->massini > 0) { 507 baryonNumberBalanceInINCL -= G4int(varntp->massini); 508 chargeNumberBalanceInINCL -= G4int(varntp->mzini); 509 // Call Fermi Break-up 510 G4double nuclearMass = G4NucleiProperties::GetNuclearMass(G4int(varntp->massini), G4int(varntp->mzini)) + varntp->exini * MeV; 511 G4LorentzVector fragmentMomentum(varntp->pxrem * MeV, varntp->pyrem * MeV, varntp->pzrem * MeV, 512 varntp->erecrem * MeV + nuclearMass); 513 G4double momentumScaling = G4InclUtils::calculate4MomentumScaling(G4int(varntp->massini), G4int(varntp->mzini), 514 varntp->exini, 515 varntp->erecrem, 516 varntp->pxrem, 517 varntp->pyrem, 518 varntp->pzrem); 519 G4LorentzVector p4(momentumScaling * varntp->pxrem * MeV, momentumScaling * varntp->pyrem * MeV, 520 momentumScaling * varntp->pzrem * MeV, 521 varntp->erecrem + nuclearMass); 522 523 // For four-momentum, baryon number and charge conservation check: 524 G4LorentzVector fourMomentumBalance = p4; 525 G4int baryonNumberBalance = G4int(varntp->massini); 526 G4int chargeBalance = G4int(varntp->mzini); 527 528 G4LorentzRotation toFragmentZ; 529 toFragmentZ.rotateZ(-p4.theta()); 530 toFragmentZ.rotateY(-p4.phi()); 531 G4LorentzRotation toFragmentLab = toFragmentZ.inverse(); 532 p4 *= toFragmentZ; 533 534 G4LorentzVector p4rest = p4; 535 p4rest.boost(-p4.boostVector()); 536 if(verboseLevel > 0) { 537 G4cout <<"Cascade remnant nucleus:" << G4endl; 538 G4cout <<"p4: " << G4endl; 539 G4cout <<" px: " << p4.px() <<" py: " << p4.py() <<" pz: " << p4.pz() << G4endl; 540 G4cout <<" E = " << p4.e() << G4endl; 541 542 G4cout <<"p4rest: " << G4endl; 543 G4cout <<" px: " << p4rest.px() <<" py: " << p4rest.py() <<" pz: " << p4rest.pz() << G4endl; 544 G4cout <<" E = " << p4rest.e() << G4endl; 545 } 546 547 G4Fragment theCascadeRemnant(G4int(varntp->massini), G4int(varntp->mzini), p4rest); 548 thePrecoResult = thePrecoModel->DeExcite(theCascadeRemnant); 549 if(thePrecoResult != 0) { 550 G4ReactionProductVector::iterator fragment; 551 for(fragment = thePrecoResult->begin(); fragment != thePrecoResult->end(); fragment++) { 552 G4ParticleDefinition *theFragmentDefinition = (*fragment)->GetDefinition(); 553 554 if(theFragmentDefinition != 0) { 555 G4DynamicParticle *theFragment = new G4DynamicParticle(theFragmentDefinition, (*fragment)->GetMomentum()); 556 G4LorentzVector labMomentum = theFragment->Get4Momentum(); 557 labMomentum.boost(p4.boostVector()); 558 labMomentum *= toFragmentLab; 559 labMomentum *= toLabFrame; 560 theFragment->Set4Momentum(labMomentum); 561 fourMomentumBalance -= theFragment->Get4Momentum(); 562 baryonNumberBalance -= theFragmentDefinition->GetAtomicMass(); 563 chargeBalance -= theFragmentDefinition->GetAtomicNumber(); 564 if(verboseLevel > 0) { 565 G4cout <<"Resulting fragment: " << G4endl; 566 G4cout <<" kinetic energy = " << theFragment->GetKineticEnergy() / MeV << " MeV" << G4endl; 567 G4cout <<" momentum = " << theFragment->GetMomentum().mag() / MeV << " MeV" << G4endl; 568 } 569 theResult.AddSecondary(theFragment); 570 } else { 571 G4cout <<"G4InclCascadeInterface: Error. Fragment produced by Fermi break-up does not exist." << G4endl; 572 G4cout <<"Resulting fragment: " << G4endl; 573 G4cout <<" momentum = " << (*fragment)->GetMomentum().mag() / MeV << " MeV" << G4endl; 574 } 575 } 576 delete thePrecoResult; 577 thePrecoResult = 0; 578 579 if(verboseLevel > 1 && std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 580 G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl; 581 G4cout <<"Magnitude: " << fourMomentumBalance.mag() / MeV << " MeV" << G4endl; 582 G4cout <<"Vector components (px, py, pz, E) = (" 583 << fourMomentumBalance.px() << ", " 584 << fourMomentumBalance.py() << ", " 585 << fourMomentumBalance.pz() << ", " 586 << fourMomentumBalance.e() << ")" << G4endl; 587 } 588 if(baryonNumberBalance != 0 && verboseLevel > 1) { 589 G4cout <<"Baryon number balance after remnant nucleus Fermi break-up: " << baryonNumberBalance << G4endl; 590 } 591 if(chargeBalance != 0 && verboseLevel > 1) { 592 G4cout <<"Charge balance after remnant nucleus Fermi break-up: " << chargeBalance << G4endl; 593 } 594 } 595 } 596 294 597 varntp->ntrack = 0; // Clean up the number of generated particles in the event. 598 599 if(baryonNumberBalanceInINCL != 0 && verboseLevel > 1) { 600 G4cout <<"Event " << eventNumber <<": G4InclLightIonInterface: Baryon number conservation problem in INCL detected!" << G4endl; 601 G4cout <<"Baryon number balance: " << baryonNumberBalanceInINCL << G4endl; 602 if(baryonNumberBalanceInINCL < 0) { 603 G4cout <<"Event " << eventNumber <<": Too many outcoming baryons!" << G4endl; 604 } else if(baryonNumberBalanceInINCL > 0) { 605 G4cout <<"Event " << eventNumber <<": Too few outcoming baryons!" << G4endl; 606 } 607 } 608 609 if(chargeNumberBalanceInINCL != 0 && verboseLevel > 1) { 610 G4cout <<"Event " << eventNumber <<": G4InclLightIonInterface: Charge number conservation problem in INCL detected!" << G4endl; 611 G4cout <<"Event " << eventNumber <<": Charge number balance: " << chargeNumberBalanceInINCL << G4endl; 612 } 295 613 } 296 614 /** … … 304 622 cascadeParticle = new G4DynamicParticle(theTableOfParticles->FindParticle(aTrack.GetDefinition()), aTrack.Get4Momentum()); 305 623 306 theResult.AddSecondary(cascadeParticle);624 result.push_back(cascadeParticle); 307 625 308 626 if(verboseLevel > 1) { … … 317 635 G4cout <<"G4InclLightIonInterface: Unknown bullet type" << G4endl; 318 636 G4cout <<"Bullet particle name: " << cascadeParticle->GetDefinition()->GetParticleName() << G4endl; 319 } 637 } 320 638 if(verboseLevel > 3) { 321 639 diagdata <<"G4InclLightIonInterface: Unknown bullet type" << G4endl; … … 324 642 } 325 643 326 if(( targetA == 1) && (targetZ== 1)) { // Unsupported target644 if((calincl->targetA() == 1) && (calincl->targetZ() == 1)) { // Unsupported target 327 645 if(verboseLevel > 1) { 328 646 G4cout <<"Unsupported target: " << G4endl; 329 G4cout <<"Target A: " << targetA<< G4endl;330 G4cout <<"TargetZ: " << targetZ<< G4endl;647 G4cout <<"Target A: " << calincl->targetA() << G4endl; 648 G4cout <<"TargetZ: " << calincl->targetZ() << G4endl; 331 649 } 332 650 if(verboseLevel > 3) { 333 651 diagdata <<"Unsupported target: " << G4endl; 334 diagdata <<"Target A: " << targetA<< G4endl;335 diagdata <<"TargetZ: " << targetZ<< G4endl;336 } 337 } 338 339 if( bulletE< 100) { // INCL does not support E < 100 MeV.652 diagdata <<"Target A: " << calincl->targetA() << G4endl; 653 diagdata <<"TargetZ: " << calincl->targetZ() << G4endl; 654 } 655 } 656 657 if(calincl->bulletE() < 100) { // INCL does not support E < 100 MeV. 340 658 if(verboseLevel > 1) { 341 G4cout <<"Unsupported bullet energy: " << bulletE<< " MeV. (Lower limit is 100 MeV)." << G4endl;659 G4cout <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl; 342 660 G4cout <<"WARNING: Returning the original bullet with original energy back to Geant4." << G4endl; 343 661 } 344 662 if(verboseLevel > 3) { 345 diagdata <<"Unsupported bullet energy: " << bulletE << " MeV. (Lower limit is 100 MeV)." << G4endl; 346 } 347 } 663 diagdata <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl; 664 } 665 } 666 348 667 if(verboseLevel > 3) { 349 668 diagdata <<"WARNING: returning the original bullet with original energy back to Geant4." << G4endl; … … 351 670 } 352 671 672 // Finally copy the accumulated secondaries into the result collection: 673 G4ThreeVector boostVector = aTrack.Get4Momentum().boostVector(); 674 G4LorentzRotation boostBack = toBreit.inverse(); 675 676 for(std::vector<G4DynamicParticle*>::iterator i = result.begin(); i != result.end(); ++i) { 677 // If the calculation was performed in inverse kinematics we have to 678 // convert the result back... 679 if(calincl->isInverseKinematics()) { 680 G4LorentzVector mom = (*i)->Get4Momentum(); 681 mom.setPz(-1.0 * mom.pz()); // Reverse the z-component of the momentum vector 682 mom *= boostBack; 683 (*i)->Set4Momentum(mom); 684 } 685 theResult.AddSecondary((*i)); 686 } 687 688 delete calincl; 689 calincl = 0; 353 690 return &theResult; 354 691 }
Note: See TracChangeset
for help on using the changeset viewer.