Ignore:
Timestamp:
Nov 25, 2009, 5:13:58 PM (15 years ago)
Author:
garnier
Message:

update CVS release candidate geant4.9.3.01

Location:
trunk/source/processes/hadronic/models/incl/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/incl/src/G4Abla.cc

    r962 r1196  
    2424// ********************************************************************
    2525//
    26 // $Id: G4Abla.cc,v 1.19 2008/09/15 08:16:45 kaitanie Exp $
     26// $Id: G4Abla.cc,v 1.20 2009/11/18 10:43:14 kaitanie Exp $
    2727// Translation of INCL4.2/ABLA V3
    2828// Pekka Kaitaniemi, HIP (translation)
    … …  
    651651      G4double bil_pz_pf1 = pf1_rem[3];
    652652      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
    653654        bil_e_pf1 = bil_e_pf1 - (std::pow(varntp->plab[ipf1],2) + std::pow(varntp->enerj[ipf1],2))/(2.0*(varntp->enerj[ipf1]));
    654655        cst = std::cos(varntp->tetlab[ipf1]/57.2957795);
    … …  
    668669      G4double bil_pz_pf2 = pf2_rem[3];
    669670      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
    670672        bil_e_pf2 = bil_e_pf2 - (std::pow(varntp->plab[ipf2],2) + std::pow(varntp->enerj[ipf2],2))/(2.0*(varntp->enerj[ipf2]));
    671673        G4double cst = std::cos(varntp->tetlab[ipf2]/57.2957795);

  • trunk/source/processes/hadronic/models/incl/src/G4Incl.cc

    r962 r1196  
    2424// ********************************************************************
    2525//
    26 // $Id: G4Incl.cc,v 1.20 2008/11/06 10:11:27 kaitanie Exp $
     26// $Id: G4Incl.cc,v 1.22 2009/11/18 10:43:14 kaitanie Exp $
    2727// Translation of INCL4.2/ABLA V3
    2828// Pekka Kaitaniemi, HIP (translation)
    … …  
    154154 *Methods for debugging.
    155155 */
     156
     157void 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
    156172G4double G4Incl::energyTest(G4int i)
    157173{
    … …  
    617633    if((std::fabs(pzbil-pbeam) > 5.0) || (std::sqrt(std::pow(pxbil,2)+std::pow(pybil,2)) >= 3.0)) {
    618634      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;
    620638      }
    621639    }
    … …  
    981999    if((std::fabs(pzbil - pbeam) > 5.0) || (std::sqrt(std::pow(pxbil,2) + std::pow(pybil,2)) >= 3.0)) {
    9821000      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;
    9841004      }
    9851005    }
    … …  
    13871407    return (saxw->y[0][saxw->imat]);
    13881408  }
    1389   else {
     1409  else { // tz > 0
    13901410    for(G4int i = 1; i < saxw->n; i++) {
    1391       j = i - 1;
     1411      j = i;
    13921412      tz = xv - saxw->x[j][saxw->imat];
    1393       if(tz < 0) {
     1413      if(tz <= 0) {
    13941414        break;
    13951415      }
    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;
    14041427}
    14051428
    … …  
    15471570  x1[3] = 793.0/720.0;
    15481571  x1[4] = 157.0/160.0;
    1549   nb = int(std::floor(((ra - ri)/step + 1.0000000001))); // 1.0000000001 -> 0.0
     1572  nb = int(std::floor(((ra - ri)/dr + 1.0000000001))); // 1.0000000001 -> 0.0
    15501573  dr = (ra - ri)/(double(nb - 1));
    15511574  res = 0.0;
    … …  
    23132336  // espace de phases test (r et p) pour pauli:
    23142337  // 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;
    23162342  G4double pbl=200.0;
    23172343
    … …  
    24522478  //G4double bimpact=b;
    24532479  bimpact = b;
    2454   G4double tnor;
     2480  G4double tnor = 0.0;
    24552481
    24562482  if(ws->nosurf != -2) { // la suite, c'est la version temps avant 2001
    … …  
    35813607  }
    35823608 
    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:
    36083624  if (bb2 < bmax2) {
    36093625    goto pnu220;
    … …  
    57405756  // backward scattering according the parametrization of ref
    57415757  // prc56(1997)1
    5742 
    5743   if(((m1+m2) != 1) || (iso == 0)) {
    5744     standardRandom(&rndm, &(hazard->igraine[7]));
    5745     apt = 1.0;
    5746     if (pl > 800.0) {
    5747       apt = std::pow((800.0/pl),2);
    5748     }
    5749     if ((iexpi == 1) || (rndm > (1./(1.+apt)))) {
    5750       ii = is1;
    5751       is1 = is2;
    5752       is2 = ii;
    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:
    57555771
    57565772  debugOutput = am(p1,p2,p3,e1);
    … …  
    61086124  }
    61096125
    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]));
    61216134  fi = 6.2832*rndm;
    61226135  cfi = std::cos(fi);
    61236136  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;
    61276140  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:
    61426153  hel = 0.0;                                                       
    61436154  w1 = q1*q1 + q2*q2 + q3*q3;
    … …  
    62506261      goto newt50;
    62516262    }
    6252     if (bl1->ta < bl5->tlg[l1]) { // tlg(12)->tlg[11]
     6263    if (bl1->ta < bl5->tlg[l2]) { // tlg(12)->tlg[11]
    62536264      goto newt50;
    62546265    }
    … …  
    65536564    xr2 = xr*xr;
    65546565    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));
    65566567    rs = std::sqrt(bl3->x1[l]*bl3->x1[l] + bl3->x2[l]*bl3->x2[l] + bl3->x3[l]*bl3->x3[l]);
    65576568    if (ws->nosurf <= 0) {
    … …  
    68726883}
    68736884
    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 
     6885G4double 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;
    68786891  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;
    68826893  G4double p = std::sqrt(t2);
    68836894  G4double r = r2;
    68846895  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.:
    68886899    xv = p/pf;
    6889     r = interpolateFunction(xv);
     6900    if(t2 <= pf2) {
     6901      r = interpolateFunction(xv);
     6902    } else {
     6903      r = ws->rmaxws;
     6904    }
    68906905    r = r*r;
    6891     if (t2 > pf2) {
    6892       r = std::pow(ws->rmaxws,2);
    6893     }
    6894   }
    6895 
     6906  }
     6907 ref21:
    68966908  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;
    69056910  t1 = x1*p1 + x2*p2 + x3*p3;
    69066911  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;
    69176927}
    69186928
    … …  
    77417751  G4int valueFloor = int(std::floor(a));
    77427752
    7743   if(std::abs(value - valueCeil) < std::abs(value - valueFloor)) {
     7753  if(std::abs(value - valueCeil) <= std::abs(value - valueFloor)) {
    77447754    return valueCeil;
    77457755  }

  • trunk/source/processes/hadronic/models/incl/src/G4InclAblaCascadeInterface.cc

    r819 r1196  
    2424// ********************************************************************
    2525//
    26 // $Id: G4InclAblaCascadeInterface.cc,v 1.10 2007/12/10 16:31:55 gunter Exp $
     26// $Id: G4InclAblaCascadeInterface.cc,v 1.12 2009/11/18 10:43:14 kaitanie Exp $
    2727// Translation of INCL4.2/ABLA V3
    2828// Pekka Kaitaniemi, HIP (translation)
    … …  
    3838#include "CLHEP/Random/Random.h"
    3939
    40 G4InclAblaCascadeInterface::G4InclAblaCascadeInterface()
     40G4InclAblaCascadeInterface::G4InclAblaCascadeInterface(const G4String& nam)
     41  :G4VIntraNuclearTransportModel(nam)
    4142{
    4243  hazard = new G4Hazard();
Note: See TracChangeset for help on using the changeset viewer.