Changeset 1196 for trunk/source/processes/hadronic/models/incl
- Timestamp:
- Nov 25, 2009, 5:13:58 PM (15 years ago)
- Location:
- trunk/source/processes/hadronic/models/incl
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/incl/History
r962 r1196 4 4 Geant4 - an Object-Oriented Toolkit for Physics Simulation 5 5 ========================================================== 6 $Id: History,v 1.1 7 2008/11/06 10:11:27kaitanie Exp $6 $Id: History,v 1.19 2009/11/18 10:43:14 kaitanie Exp $ 7 7 --------------------------------------------------------------------- 8 8 … … 16 16 * Please list in reverse chronological order (last date on top) 17 17 --------------------------------------------------------------- 18 19 18 November 2009 - Pekka Kaitaniemi (hadr-incl-V09-02-01) 20 --------------------------------------------------------- 21 - Added safeguard division by zero (or negative) energy in ABLA 22 fission fragment handling 23 - Minor fix in the INCL particle reflection time calculation 24 25 17 October 2009 - Pekka Kaitaniemi (hadr-incl-V09-02-00) 26 -------------------------------------------------------- 27 - Bugfixes to INCL nuclear potential handling: 28 The interpolation function produced discontinuous results. This was 29 due to an array index off-by-one bug. This bugfix corrects the 30 impact parameter distribution. 31 - Fixed several variable initialization issues in INCL. 32 - Fixed severalFORTRAN to C++ translation issues in INCL. 18 33 19 34 06 November 2008 - Pekka Kaitaniemi (hadr-incl-V09-01-05) -
trunk/source/processes/hadronic/models/incl/include/G4Incl.hh
r962 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4Incl.hh,v 1.1 3 2008/06/25 17:20:04 kaitanie Exp $26 // $Id: G4Incl.hh,v 1.15 2009/11/18 10:43:14 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 78 78 ~G4Incl(); // Destructor 79 79 80 void dumpParticles(); 80 81 G4double energyTest(G4int i); // Test for NaN energy of particle i. 81 82 void dumpBl5(std::ofstream& dumpOut); // Dump the contents of G4Bl5. … … 474 475 * @return a double value 475 476 */ 476 G4double ref(G4double x1, G4double x2, G4doublex3, G4double p1,477 G4double ref(G4double &x1, G4double &x2, G4double &x3, G4double p1, 477 478 G4double p2, G4double p3, G4double E, G4double r2); 478 479 -
trunk/source/processes/hadronic/models/incl/include/G4InclAblaCascadeInterface.hh
r819 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclAblaCascadeInterface.hh,v 1. 6 2007/10/31 10:44:22 miheikkiExp $26 // $Id: G4InclAblaCascadeInterface.hh,v 1.8 2009/11/18 10:43:14 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 66 66 67 67 /** 68 * <h1>INCL intra-nuclear cascade with built-in ABLA de-excitation</h1> 69 * 68 70 * Interface for INCL/ABLA. This interface handles basic hadron 69 71 * bullet particles (protons, neutrons, pions). 72 * 73 * Example usage in case of protons: 74 * @code 75 * G4InclAblaCascadeInterface* inclModel = new G4InclAblaCascadeInterface; 76 * inclModel -> SetMinEnergy(0.0 * MeV); // Set the energy limits 77 * inclModel -> SetMaxEnergy(3.0 * GeV); 78 * 79 * G4ProtonInelasticProcess* protonInelasticProcess = new G4ProtonInelasticProcess(); 80 * G4ProtonInelasticCrossSection* protonInelasticCrossSection = new G4ProtonInelasticCrossSection(); 81 * 82 * protonInelasticProcess -> RegisterMe(inclModel); 83 * protonInelasticProcess -> AddDataSet(protonInelasticCrossSection); 84 * 85 * particle = G4Proton::Proton(); 86 * processManager = particle -> GetProcessManager(); 87 * processManager -> AddDiscreteProcess(protonInelasticProcess); 88 * @endcode 89 * The same setup procedure is needed for neutron and pion inelastic processes 90 * as well. 91 * 70 92 * @see G4InclAblaLightIonInterface 71 93 */ … … 77 99 * Basic constructor. 78 100 */ 79 G4InclAblaCascadeInterface( );101 G4InclAblaCascadeInterface(const G4String& name = "INCL/ABLA Cascade"); 80 102 81 103 -
trunk/source/processes/hadronic/models/incl/include/G4InclDataDefs.hh
r962 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclDataDefs.hh,v 1. 5 2008/06/25 17:20:04 kaitanie Exp $26 // $Id: G4InclDataDefs.hh,v 1.7 2009/11/18 10:43:14 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 323 323 ~G4Bl2() {}; 324 324 325 void dump() { 326 G4cout <<"Avatars: (number of avatars = " << k << ")" << G4endl; 327 for(G4int i = 0; i <= k; i++) { 328 G4cout <<"i = " << i << G4endl; 329 G4cout <<"crois[" << i << "] = " << crois[i] << G4endl; 330 G4cout <<"ind[" << i << "] = " << ind[i] << G4endl; 331 G4cout <<"jnd[" << i << "] = " << jnd[i] << G4endl; 332 } 333 } 334 325 335 /** 326 336 * -
trunk/source/processes/hadronic/models/incl/src/G4Abla.cc
r962 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4Abla.cc,v 1. 19 2008/09/15 08:16:45kaitanie Exp $26 // $Id: G4Abla.cc,v 1.20 2009/11/18 10:43:14 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 651 651 G4double bil_pz_pf1 = pf1_rem[3]; 652 652 for(G4int ipf1 = lmi_pf1; ipf1 <= lma_pf1; ipf1++) { //do ipf1=lmi_pf1,lma_pf1 653 if(varntp->enerj[ipf1] <= 0.0) continue; // Safeguard against a division by zero 653 654 bil_e_pf1 = bil_e_pf1 - (std::pow(varntp->plab[ipf1],2) + std::pow(varntp->enerj[ipf1],2))/(2.0*(varntp->enerj[ipf1])); 654 655 cst = std::cos(varntp->tetlab[ipf1]/57.2957795); … … 668 669 G4double bil_pz_pf2 = pf2_rem[3]; 669 670 for(G4int ipf2 = lmi_pf2; ipf2 <= lma_pf2; ipf2++) { //do ipf2=lmi_pf2,lma_pf2 671 if(varntp->enerj[ipf2] <= 0.0) continue; // Safeguard against a division by zero 670 672 bil_e_pf2 = bil_e_pf2 - (std::pow(varntp->plab[ipf2],2) + std::pow(varntp->enerj[ipf2],2))/(2.0*(varntp->enerj[ipf2])); 671 673 G4double cst = std::cos(varntp->tetlab[ipf2]/57.2957795); -
trunk/source/processes/hadronic/models/incl/src/G4Incl.cc
r962 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4Incl.cc,v 1.2 0 2008/11/06 10:11:27kaitanie Exp $26 // $Id: G4Incl.cc,v 1.22 2009/11/18 10:43:14 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 154 154 *Methods for debugging. 155 155 */ 156 157 void G4Incl::dumpParticles() 158 { 159 G4int ia = bl3->ia1 + bl3->ia2; 160 G4cout <<"Nucleons: (number of nucleons = " << ia << ")" << G4endl; 161 for(G4int i = 0; i <= ia; i++) { 162 G4cout <<"x1(" << i << ") = " << bl3->x1[i] << G4endl; 163 G4cout <<"x2(" << i << ") = " << bl3->x2[i] << G4endl; 164 G4cout <<"x3(" << i << ") = " << bl3->x3[i] << G4endl; 165 G4cout <<"p1(" << i << ") = " << bl1->p1[i] << G4endl; 166 G4cout <<"p2(" << i << ") = " << bl1->p2[i] << G4endl; 167 G4cout <<"p3(" << i << ") = " << bl1->p3[i] << G4endl; 168 G4cout <<"eps(" << i << ") = " << bl1->eps[i] << G4endl; 169 } 170 } 171 156 172 G4double G4Incl::energyTest(G4int i) 157 173 { … … 617 633 if((std::fabs(pzbil-pbeam) > 5.0) || (std::sqrt(std::pow(pxbil,2)+std::pow(pybil,2)) >= 3.0)) { 618 634 if(verboseLevel > 3) { 619 G4cout <<"bad momentum conservation after incl:" << G4endl; 635 G4cout <<"Bad momentum conservation after INCL:" << G4endl; 636 G4cout <<"delta Pz = " << std::fabs(pzbil - pbeam) << G4endl; 637 G4cout <<" Pt = " << std::sqrt(std::pow(pxbil, 2) + std::pow(pybil, 2)) << G4endl; 620 638 } 621 639 } … … 981 999 if((std::fabs(pzbil - pbeam) > 5.0) || (std::sqrt(std::pow(pxbil,2) + std::pow(pybil,2)) >= 3.0)) { 982 1000 if(verboseLevel > 3) { 983 G4cout <<"bad momentum conservation after incl:" << G4endl; 1001 G4cout <<"Bad momentum conservation after INCL:" << G4endl; 1002 G4cout <<"delta Pz = " << std::fabs(pzbil - pbeam) << G4endl; 1003 G4cout <<" Pt = " << std::sqrt(std::pow(pxbil, 2) + std::pow(pybil, 2)) << G4endl; 984 1004 } 985 1005 } … … 1387 1407 return (saxw->y[0][saxw->imat]); 1388 1408 } 1389 else { 1409 else { // tz > 0 1390 1410 for(G4int i = 1; i < saxw->n; i++) { 1391 j = i - 1;1411 j = i; 1392 1412 tz = xv - saxw->x[j][saxw->imat]; 1393 if(tz < 0) {1413 if(tz <= 0) { 1394 1414 break; 1395 1415 } 1396 else if(tz == 0) { 1397 return saxw->y[j][saxw->imat]; 1398 } 1399 } 1400 } 1401 1402 G4double dgx = xv - saxw->x[j][saxw->imat]; 1403 return(saxw->y[j][saxw->imat] + saxw->s[j][saxw->imat]*dgx); 1416 } 1417 if(tz >= 0) { 1418 return saxw->y[j][saxw->imat]; 1419 } else if(tz < 0.0) { 1420 j = j - 1; 1421 G4double dgx = xv - saxw->x[j][saxw->imat]; 1422 return(saxw->y[j][saxw->imat] + saxw->s[j][saxw->imat]*dgx); 1423 } 1424 } 1425 1426 return 0.0; 1404 1427 } 1405 1428 … … 1547 1570 x1[3] = 793.0/720.0; 1548 1571 x1[4] = 157.0/160.0; 1549 nb = int(std::floor(((ra - ri)/ step+ 1.0000000001))); // 1.0000000001 -> 0.01572 nb = int(std::floor(((ra - ri)/dr + 1.0000000001))); // 1.0000000001 -> 0.0 1550 1573 dr = (ra - ri)/(double(nb - 1)); 1551 1574 res = 0.0; … … 2313 2336 // espace de phases test (r et p) pour pauli: 2314 2337 // valeur recommandee par j.c. v-test=0.589 h**3: 2315 G4double rbl = 2.0; 2338 // G4double rbl = 2.0; 2339 2340 // Valeur pour avoir V-test=2.38 h**3 (avec pbl=200) 2341 G4double rbl=3.1848; 2316 2342 G4double pbl=200.0; 2317 2343 … … 2452 2478 //G4double bimpact=b; 2453 2479 bimpact = b; 2454 G4double tnor ;2480 G4double tnor = 0.0; 2455 2481 2456 2482 if(ws->nosurf != -2) { // la suite, c'est la version temps avant 2001 … … 3581 3607 } 3582 3608 3583 // Replaced goto structure: 3584 // if (k3 == 1) go to 260 3585 // if (k4 == 0) go to 260 3586 if(k3 != 1 && k4 != 0) { 3587 mg=bl1->ind1[bl9->l1]+bl1->ind1[bl9->l2]; 3588 isos=bl1->ind2[bl9->l1]+bl1->ind2[bl9->l2]; 3589 } 3590 if((k3 != 1) && (k4 != 0) && (mg == 1)) { 3591 // if (mg != 1) go to 260 3592 ldel = bl9->l2; 3593 if(mg-bl1->ind1[bl9->l1] == 0) { 3594 ldel = bl9->l1; 3595 } 3596 bl6->xx10 = std::sqrt(std::pow(bl1->eps[ldel],2) - std::pow(bl1->p1[ldel],2) - std::pow(bl1->p2[ldel],2) - std::pow(bl1->p3[ldel],2)); 3597 bl6->isa = bl1->ind2[ldel]; 3598 bmax2 = totalCrossSection(sq,mg,isos)/31.415926; 3599 if (k5 == 0 && mg != 0) { 3600 bmax2 = bmax2 - lowEnergy(sq,mg,isos)/31.415926; 3601 } 3602 // go to 261 3603 } 3604 else { 3605 bmax2 = totalCrossSection(sq,mg,isos)/31.41592; 3606 } 3607 3609 if (k3 == 1) goto pnu260; 3610 if (k4 == 0) goto pnu260; 3611 mg = bl1->ind1[bl9->l1] + bl1->ind1[bl9->l2]; 3612 isos = bl1->ind2[bl9->l1] + bl1->ind2[bl9->l2]; 3613 if (mg != 1) goto pnu260; 3614 ldel = bl9->l2; 3615 if(mg - bl1->ind1[bl9->l1] == 0) ldel = bl9->l1; 3616 bl6->xx10 = std::sqrt(std::pow(bl1->eps[ldel],2) - std::pow(bl1->p1[ldel], 2) - std::pow(bl1->p2[ldel], 2) - std::pow(bl1->p3[ldel], 2)); 3617 bl6->isa = bl1->ind2[ldel]; 3618 bmax2 = totalCrossSection(sq,mg,isos)/31.415926; 3619 if (k5 == 0 && mg != 0) bmax2 = bmax2 - lowEnergy(sq,mg,isos)/31.415926; 3620 goto pnu261; 3621 pnu260: 3622 bmax2 = totalCrossSection(sq,mg,isos)/31.41592; 3623 pnu261: 3608 3624 if (bb2 < bmax2) { 3609 3625 goto pnu220; … … 5740 5756 // backward scattering according the parametrization of ref 5741 5757 // prc56(1997)1 5742 5743 if (((m1+m2) != 1) || (iso == 0)) {5744 standardRandom(&rndm,&(hazard->igraine[7]));5745 5746 5747 apt = std::pow((800.0/pl),2);5748 }5749 if ((iexpi == 1) || (rndm > (1./(1.+apt)))) {5750 5751 5752 5753 }5754 }5758 if (m1+m2 == 1) goto collis133; 5759 if (iso != 0) goto collis133; 5760 standardRandom(&rndm,&(hazard->igraine[7])); 5761 apt = 1.0; 5762 if (pl > 800.0) { 5763 apt = std::pow(800.0/pl,2); 5764 } //endif 5765 if (iexpi == 1 || rndm > 1.0/(1.0+apt)) { // then 5766 ii = is1; 5767 is1 = is2; 5768 is2 = ii; 5769 } // endif 5770 collis133: 5755 5771 5756 5772 debugOutput = am(p1,p2,p3,e1); … … 6108 6124 } 6109 6125 6110 do { 6111 standardRandom(&rndm, &(hazard->igraine[7])); 6112 ctet = -1.0 + 2.0*rndm; 6113 if(std::fabs(ctet) > 1.0) { 6114 ctet = sign(1.0,ctet); 6115 } 6116 stet = std::sqrt(1.0 - std::pow(ctet,2)); 6117 standardRandom(&rndm, &(hazard->igraine[9])); 6118 } while(rndm > ((1.0 + 3.0*hel*std::pow(ctet,2))/(1.0 + 3.0*hel))); 6119 6120 standardRandom(&rndm, &(hazard->igraine[8])); 6126 decay2100: 6127 standardRandom(&rndm,&(hazard->igraine[7])); 6128 ctet = -1.0 + 2.0*rndm; 6129 if(std::abs(ctet) > 1.0) ctet = sign(1.0,ctet); 6130 stet = std::sqrt(1.0 - std::pow(ctet, 2)); 6131 standardRandom(&rndm,&(hazard->igraine[9])); 6132 if (rndm > ((1.0 + 3.0 * hel * std::pow(ctet,2))/(1.0 + 3.0*hel))) goto decay2100; 6133 standardRandom(&rndm,&(hazard->igraine[8])); 6121 6134 fi = 6.2832*rndm; 6122 6135 cfi = std::cos(fi); 6123 6136 sfi = std::sin(fi); 6124 beta = std::sqrt(b1*b1 +b2*b2+b3*b3);6125 6126 sal = std::sqrt(std::pow(b1, 2) + std::pow(b2,2))/beta;6137 beta = std::sqrt(b1*b1 + b2*b2 + b3*b3); 6138 if (beta < 1.0e-10) goto decay2101; 6139 sal = std::sqrt(std::pow(b1, 2) + std::pow(b2, 2))/beta; 6127 6140 cal = b3/beta; 6128 6129 if((beta >= 1.0e-10) || (sal >= 1.0e-6)) { 6130 t1 = ctet + cal*stet*sfi/sal; 6131 t2 = stet/sal; 6132 q1 = xq*(b1*t1 + b2*t2*cfi)/beta; 6133 q2 = xq*(b2*t1 - b1*t2*cfi)/beta; 6134 q3 = xq*(b3*t1/beta - t2*sfi); 6135 } 6136 else { 6137 q1 = xq*stet*cfi; 6138 q2 = xq*stet*sfi; 6139 q3 = xq*ctet; 6140 } 6141 6141 if (sal < 1.0e-6) goto decay2101; 6142 t1 = ctet + cal*stet*sfi/sal; 6143 t2 = stet/sal; 6144 q1 = xq*(b1*t1 + b2*t2*cfi)/beta; 6145 q2 = xq*(b2*t1 - b1*t2*cfi)/beta; 6146 q3 = xq*(b3*t1/beta - t2*sfi); 6147 goto decay2102; 6148 decay2101: 6149 q1 = xq * stet*cfi; 6150 q2 = xq * stet*sfi; 6151 q3 = xq * ctet; 6152 decay2102: 6142 6153 hel = 0.0; 6143 6154 w1 = q1*q1 + q2*q2 + q3*q3; … … 6250 6261 goto newt50; 6251 6262 } 6252 if (bl1->ta < bl5->tlg[l 1]) { // tlg(12)->tlg[11]6263 if (bl1->ta < bl5->tlg[l2]) { // tlg(12)->tlg[11] 6253 6264 goto newt50; 6254 6265 } … … 6553 6564 xr2 = xr*xr; 6554 6565 pr2 = pr*pr; 6555 vol = std::pow((40.0*3.1415926/3.0),2) * (std::pow((xr*pr)/(2.0*3.1415926*197. 13),3));6566 vol = std::pow((40.0*3.1415926/3.0),2) * (std::pow((xr*pr)/(2.0*3.1415926*197.33),3)); 6556 6567 rs = std::sqrt(bl3->x1[l]*bl3->x1[l] + bl3->x2[l]*bl3->x2[l] + bl3->x3[l]*bl3->x3[l]); 6557 6568 if (ws->nosurf <= 0) { … … 6872 6883 } 6873 6884 6874 G4double G4Incl::ref(G4double x1, G4double x2, G4double x3, G4double p1, G4double p2, G4double p3, G4double E, G4double r2) 6875 { 6876 const G4double pf = 270.339 , pf2 = 73083.4; 6877 6885 G4double G4Incl::ref(G4double &x1, G4double &x2, G4double &x3, G4double p1, G4double p2, G4double p3, G4double E, G4double r2) 6886 { 6887 // Surface : modif de REF 6888 // REF=TIME NECESSARY FOR A NUCLEON TO REACH THE SURFACE 6889 6890 const G4double pf = 270.33936, pf2 = 73083.4; 6878 6891 G4double ref = 0.0; 6879 G4double t1 = 0.0, t3 = 0.0, t4 = 0.0, t5 = 0.0; 6880 6881 G4double t2 = p1*p1 + p2*p2 + p3*p3; 6892 G4double t2 = p1*p1 +p2*p2 + p3*p3; 6882 6893 G4double p = std::sqrt(t2); 6883 6894 G4double r = r2; 6884 6895 G4double xv = 0.0; 6885 G4double s = 0.0;6886 6887 if (ws->nosurf <= 0) { 6896 G4double s_l = 0.0; 6897 G4double t1 = 0.0, t3 = 0.0, t4 = 0.0, t5 = 0.0; 6898 if (ws->nosurf <= 0) { // modif pour w.s.: 6888 6899 xv = p/pf; 6889 r = interpolateFunction(xv); 6900 if(t2 <= pf2) { 6901 r = interpolateFunction(xv); 6902 } else { 6903 r = ws->rmaxws; 6904 } 6890 6905 r = r*r; 6891 if (t2 > pf2) { 6892 r = std::pow(ws->rmaxws,2); 6893 } 6894 } 6895 6906 } 6907 ref21: 6896 6908 t4 = x1*x1 + x2*x2 + x3*x3; 6897 while(t4 > r) { 6898 s = std::sqrt(r*0.99/t4); 6899 x1 = x1*s; 6900 x2 = x2*s; 6901 x3 = x3*s; 6902 t4 = x1*x1 + x2*x2 + x3*x3; 6903 } 6904 6909 if (t4 > r) goto ref2; 6905 6910 t1 = x1*p1 + x2*p2 + x3*p3; 6906 6911 t3 = t1/t2; 6907 6908 t5 = t3*t3 + (r - t4)/t2; 6909 if (t5 > 0) { 6910 ref = (-t3 + std::sqrt(t5))*E; 6911 return ref; 6912 } 6913 else { 6914 ref = 10000.0; 6915 return ref; 6916 } 6912 t5 = t3*t3 + (r-t4)/t2; 6913 if (t5 > 0) goto ref1; 6914 ref = 10000.0; 6915 return ref; 6916 ref1: 6917 ref = (-t3 + std::sqrt(t5))*E; 6918 return ref; 6919 ref2: 6920 s_l = std::sqrt(r*0.99/t4); 6921 x1 = x1*s_l; 6922 x2 = x2*s_l; 6923 x3 = x3*s_l; 6924 goto ref21; 6925 6926 return 0.0; 6917 6927 } 6918 6928 … … 7741 7751 G4int valueFloor = int(std::floor(a)); 7742 7752 7743 if(std::abs(value - valueCeil) < std::abs(value - valueFloor)) {7753 if(std::abs(value - valueCeil) <= std::abs(value - valueFloor)) { 7744 7754 return valueCeil; 7745 7755 } -
trunk/source/processes/hadronic/models/incl/src/G4InclAblaCascadeInterface.cc
r819 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclAblaCascadeInterface.cc,v 1.1 0 2007/12/10 16:31:55 gunterExp $26 // $Id: G4InclAblaCascadeInterface.cc,v 1.12 2009/11/18 10:43:14 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 38 38 #include "CLHEP/Random/Random.h" 39 39 40 G4InclAblaCascadeInterface::G4InclAblaCascadeInterface() 40 G4InclAblaCascadeInterface::G4InclAblaCascadeInterface(const G4String& nam) 41 :G4VIntraNuclearTransportModel(nam) 41 42 { 42 43 hazard = new G4Hazard();
Note: See TracChangeset
for help on using the changeset viewer.