Changeset 962 for trunk/source/processes/hadronic/models/incl/src

Timestamp:

Apr 6, 2009, 12:30:29 PM (15 years ago)

Author:

garnier

Message:

update processes

Location:

trunk/source/processes/hadronic/models/incl/src

Files:

• G4Abla.cc (modified) (129 diffs)
• G4AblaEvaporation.cc (modified) (3 diffs)
• G4Incl.cc (modified) (195 diffs)
• G4InclAblaDataFile.cc (modified) (3 diffs)

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4Abla.cc,v 1.14 2007/12/03 19:36:06 miheikki Exp $ 
+// $Id: G4Abla.cc,v 1.19 2008/09/15 08:16:45 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -36 +36 @@
 #include "G4InclAblaDataFile.hh"
 #include "Randomize.hh"
-#include <assert.h>
+#include "G4InclRandomNumbers.hh"
+#include "G4Ranecu.hh"
+#include "G4AblaFissionSimfis18.hh"
+#include "G4AblaFission.hh"
 
 G4Abla::G4Abla()
 {
-
+  verboseLevel = 0;
+  ilast = 0;
 }
 
@@ -46 +50 @@
 {
   verboseLevel = 0;
-
+  ilast = 0;
   volant = volant; // ABLA internal particle data
   volant->iv = 0;
   hazard = hazard; // Random seeds
 
+  randomGenerator = new G4InclGeant4Random();
+  //randomGenerator = new G4Ranecu();
   varntp = new G4VarNtp();
   pace = new G4Pace();
   ald = new G4Ald();
-  ablamain = new G4Ablamain();
-  emdpar = new G4Emdpar();
   eenuc = new G4Eenuc();
   ec2sub = new G4Ec2sub();
@@ -67 +71 @@
 {
   verboseLevel = 0;
-  
+  ilast = 0;
   volant = aVolant; // ABLA internal particle data
   volant->iv = 0;
@@ -74 +78 @@
   varntp->ntrack = 0;
  
+  randomGenerator = new G4InclGeant4Random();
+  //randomGenerator = new G4Ranecu();
+
+  // ABLA fission
+  //  fissionModel = new G4AblaFissionSimfis18(hazard, randomGenerator);
+  fissionModel = new G4AblaFission(hazard, randomGenerator);
+  if(verboseLevel > 0) {
+    fissionModel->about();
+  }
+  verboseLevel = 0;
   pace = new G4Pace();
   ald = new G4Ald();
-  ablamain = new G4Ablamain();
-  emdpar = new G4Emdpar();
   eenuc = new G4Eenuc();
   ec2sub = new G4Ec2sub();
@@ -88 +100 @@
 G4Abla::~G4Abla()
 {
+  delete randomGenerator;
   delete pace;
   delete ald;
-  delete ablamain;
-  delete emdpar;
   delete eenuc;
   delete ec2sub;
@@ -111 +122 @@
                        G4int eventnumber)
 {
-  // C *************************  EVAPORATION KHS *********       
-
   const G4double uma = 931.4942;
   const G4double melec = 0.511;
@@ -118 +127 @@
   const G4double fmn = 939.56563;
 
-  // Rotation matrix...
-  G4double sitet = 0.0;
-  G4double R[4][4];
-
-  G4double plab1, gam1, eta1;
-
-  G4double plab2, gam2, eta2;
-
-  G4double stet1;
-  G4double stet2;
-  
-  G4int nbpevap;
-  G4int mempaw = 0, memiv = 0;
-
+  G4double alrem = 0.0, berem = 0.0, garem = 0.0;
+  G4double R[4][4]; // Rotation matrix
   G4double csdir1[4];
   G4double csdir2[4];
-
   G4double csrem[4];
-  G4double alrem = 0.0, berem = 0.0, garem = 0.0;
-
   G4double pfis_rem[4];
   G4double pf1_rem[4];
+  for(G4int init_i = 0; init_i < 4; init_i++) {
+    csdir1[init_i] = 0.0;
+    csdir2[init_i] = 0.0;
+    csrem[init_i] = 0.0;
+    pfis_rem[init_i] = 0.0;
+    pf1_rem[init_i] = 0.0;
+    for(G4int init_j = 0; init_j < 4; init_j++) {
+      R[init_i][init_j] = 0.0;
+    }
+  }
+
+  G4double plab1 = 0.0, gam1 = 0.0, eta1 = 0.0;
+  G4double plab2 = 0.0, gam2 = 0.0, eta2 = 0.0;
+
+  G4double sitet = 0.0;
+  G4double stet1 = 0.0;
+  G4double stet2 = 0.0;
+  
+  G4int nbpevap = 0;
+  G4int mempaw = 0, memiv = 0;
 
   G4double e_evapo = 0.0;
-  G4double el;
-  G4double fmcv;
-
-  G4double aff1;
-  G4double zff1;
-  G4double eff1;
-    
-  G4double aff2;
-  G4double zff2;
-  G4double eff2;
-
-  G4double v1, v2;
+  G4double el = 0.0;
+  G4double fmcv = 0.0;
+
+  G4double aff1 = 0.0;
+  G4double zff1 = 0.0;
+  G4double eff1 = 0.0;
+  G4double aff2 = 0.0;
+  G4double zff2 = 0.0;
+  G4double eff2 = 0.0;
+
+  G4double v1 = 0.0, v2 = 0.0;
 
   G4double t2 = 0.0;
-  G4double ctet1;
+  G4double ctet1 = 0.0;
   G4double ctet2 = 0.0;
-  G4double phi1;
+  G4double phi1 = 0.0;
   G4double phi2 = 0.0;
   G4double p2 = 0.0;
-  G4double epf2_out = 0.0 ; ///AH adding initialization
+  G4double epf2_out = 0.0 ;
   G4int lma_pf1 = 0, lmi_pf1 = 0;
   G4int lma_pf2 = 0, lmi_pf2 = 0;
   G4int nopart = 0;
 
-  G4double cst, sst, csf, ssf;
+  G4double cst = 0.0, sst = 0.0, csf = 0.0, ssf = 0.0;
   
-  G4double zf, af, mtota, pleva, pxeva, pyeva;
+  G4double zf = 0.0, af = 0.0, mtota = 0.0, pleva = 0.0, pxeva = 0.0, pyeva = 0.0;
   G4int ff = 0;
   G4int inum = eventnumber;
-  G4int inttype;
+  G4int inttype = 0;
   G4double esrem = excitationEnergy;
   
@@ -181 +194 @@
 
   G4double erecrem = recoilEnergy;
-  G4double trem;
+  G4double trem = 0.0;
   G4double pxrem = momX;
   G4double pyrem = momY;
   G4double pzrem = momZ;
 
-  G4double remmass;
+  G4double remmass = 0.0;
   
+  volant->clear(); // Clean up an initialize ABLA output.
+  varntp->clear(); // Clean up an initialize ABLA output.
   varntp->ntrack = 0;
-  volant->iv = 0;
+    volant->iv = 0;
+  //volant->iv = 1;
   
   G4double pcorem = std::sqrt(erecrem*(erecrem +2.*938.2796*nucleusA));
     // G4double pcorem = std::sqrt(std::pow(momX,2) + std::pow(momY,2) + std::pow(momZ,2));
-    // assert(isnan(pcorem) == false);
-  if(esrem >= 1.0e-3) { //then         
-    //   void evapora(G4double zprf, G4double aprf, G4double ee, G4double jprf, 
-    //         G4double *zf_par, G4double *af_par, G4double *mtota_par,
-    //         G4double *pleva_par, G4double *pxeva_par, G4double *pyeva_par,
-    //         G4double *ff_par, G4int *inttype_par, G4int *inum_par);
-    //     G4cout <<"Evaporating nucleus: " << G4endl;
-    //     G4cout <<"A = " << aprf << " Z = " << zprf << G4endl;
+  if(pcorem != 0) { // Guard against division by zero.
+    alrem = pxrem/pcorem;
+    berem = pyrem/pcorem;
+    garem = pzrem/pcorem;
+  } else {
+    alrem = 0.0;
+    berem = 0.0;
+    garem = 0.0;
+  }
+
+  G4int idebug = 0;
+  G4int itest = 0;
+  if(idebug == 1) {
+    zprf =   81.;
+    aprf =   201.;
+    //    ee =   86.5877686;
+    ee = 300.0;
+    jprf =   32.;
+    zf =   0.;
+    af =   0.;
+    mtota =   0.;
+    pleva =   0.;
+    pxeva =   0.;
+    pyeva =   0.;
+    ff =  -1;
+    inttype =  0;
+    inum =  2;
+  }
+  if(itest == 1) {
+    G4cout <<" PK::: EVAPORA event  " << inum << G4endl;
+    G4cout <<" PK:::    zprf =   " << zprf << G4endl;
+    G4cout <<" PK:::    aprf =   " << aprf << G4endl;
+    G4cout <<" PK:::      ee =   " << ee << G4endl;
+    G4cout <<" PK:::      eeDiff =   " << ee - 86.5877686 << G4endl;
+    G4cout <<" PK:::    jprf =   " << jprf << G4endl;
+    G4cout <<" PK:::      zf =   " << zf << G4endl;
+    G4cout <<" PK:::      af =   " << af << G4endl;
+    G4cout <<" PK:::   mtota =   " << mtota << G4endl;
+    G4cout <<" PK:::   pleva =   " << pleva << G4endl;
+    G4cout <<" PK:::   pxeva =   " << pxeva << G4endl;
+    G4cout <<" PK:::   pyeva =   " << pyeva << G4endl;
+    G4cout <<" PK:::      ff =  " << ff << G4endl;
+    G4cout <<" PK::: inttype =  " << inttype << G4endl;
+    G4cout <<" PK:::    inum =  " << inum << G4endl;
+  }
+
+  if(esrem >= 1.0e-3) {
     evapora(zprf,aprf,ee,jprf, &zf, &af, &mtota, &pleva, &pxeva, &pyeva, &ff, &inttype, &inum);
-    // assert(isnan(pleva) == false);
-    // assert(isnan(pxeva) == false);
-    // assert(isnan(pyeva) == false);
   }
   else {
@@ -213 +265 @@
     pyeva = pyrem;
     pleva = pzrem;
-  } //   endif
-  // assert(isnan(zf) == false);
-  // assert(isnan(af) == false);
-  // assert(isnan(ee) == false);
-  //                                                                       
-  // AFP,ZFP is the final fragment if no fission occurs (FF=0)                   
-  // In case of fission (FF=1) it is the nucleus that undergoes fission.          
-//   G4double zfp = idnint(zf);                                                  
-//   G4double afp = idnint(af);
-
-  if (ff == 1) { //then   
-    // ---------------------  Here, a FISSION occures --------------------------
-    //                                                                        
-    //  FEE: (EE) energy of fissioning nucleus ABOVE the fission barrier.          
-    // 
-    //   calcul des impulsions des particules evaporees (avant fission) 
-    //                 dans le systeme labo:
-
+  }
+
+  if (ff == 1) {
+    // Fission:
+    // variable ee: Energy of fissioning nucleus above the fission barrier.          
+    // Calcul des impulsions des particules evaporees (avant fission) 
+    // dans le systeme labo.
+
+    if(verboseLevel > 2)
+      G4cout << __FILE__ << ":" << __LINE__ << " Entering fission code." << G4endl;
     trem = double(erecrem);
     remmass = pace2(aprf,zprf) + aprf*uma - zprf*melec; // canonic
-    remmass = mcorem  + double(esrem);                  // ok
-    remmass = mcorem;                                   //cugnon
+//     remmass = mcorem  + double(esrem);                       // ok
+//     remmass = mcorem;                                        //cugnon
     varntp->kfis = 1;
     G4double gamrem = (remmass + trem)/remmass;
     G4double etrem = std::sqrt(trem*(trem + 2.0*remmass))/remmass;
-    // assert(isnan(etrem) == false);
+
     //  This is not treated as accurately as for the non fission case for which
     //  the remnant mass is computed to satisfy the energy conservation 
@@ -247 +291 @@
     mglms(aprf,zprf,0,&el);
     remmass = zprf*fmp + (aprf-zprf)*fmn + el + double(esrem);
-
     gamrem = std::sqrt(std::pow(pcorem,2) + std::pow(remmass,2))/remmass;
-    // assert(isnan(gamrem) == false);
     etrem = pcorem/remmass;
-    // assert(isnan(etrem) == false);
     
-    alrem = pxrem/pcorem;
-    // assert(isnan(alrem) == false);
-    berem = pyrem/pcorem;
-    // assert(isnan(berem) == false);
-    garem = pzrem/pcorem;
-    // assert(isnan(garem) == false);
-
     csrem[0] = 0.0; // Should not be used.
     csrem[1] = alrem;
     csrem[2] = berem;
     csrem[3] = garem;
-      
+    if(verboseLevel > 1) {
+      G4cout << __FILE__ << ":" << __LINE__ << " momRem = (" << pxrem << ", " << pyrem << ", " << pzrem << ")" << G4endl;
+      G4cout << __FILE__ << ":" << __LINE__ << " csrem = (" << csrem[1] << ", " << csrem[2] << ", " << csrem[3] << ")" << G4endl;
+    }
+
     // C Pour Vérif Remnant = evapo(Pre fission) + Noyau_fissionant (systÚme  Remnant)
     G4double bil_e = 0.0;
@@ -273 +311 @@
     
     for(G4int iloc = 1; iloc <= volant->iv; iloc++) { //DO iloc=1,iv
-      // assert(isnan(volant->zpcv[iloc]) == false);
-      //      assert(volant->acv[iloc] != 0);
-      //      assert(volant->zpcv[iloc] != 0);
       mglms(double(volant->acv[iloc]),double(volant->zpcv[iloc]),0,&el);
-      // assert(isnan(el) == false);
       masse = volant->zpcv[iloc]*fmp + (volant->acv[iloc] - volant->zpcv[iloc])*fmn + el;
-      // assert(isnan(masse) == false);
       bil_e = bil_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2));
-      // assert(isnan(bil_e) == false); 
       bil_px = bil_px + volant->pcv[iloc]*(volant->xcv[iloc]);
       bil_py = bil_py + volant->pcv[iloc]*(volant->ycv[iloc]);
@@ -315 +347 @@
     //                       G4double &a2,G4double &z2,G4double &e2,G4double &v2);
 
-    fissionDistri(af,zf,ee,aff1,zff1,eff1,v1,aff2,zff2,eff2,v2);
-
+    //fissionModel->fissionDistri(af,zf,ee,aff1,zff1,eff1,v1,aff2,zff2,eff2,v2);
+    fissionModel->doFission(af,zf,ee,aff1,zff1,eff1,v1,aff2,zff2,eff2,v2);
     // C verif des A et Z decimaux:
     G4int na_f = int(std::floor(af + 0.5));
@@ -342 +374 @@
     //G4double ef = fb->efa[kboud][jboud]; // barriere de fission
     G4double ef = fb->efa[jboud][kboud]; // barriere de fission
-    // assert(isnan(ef) == false);
     varntp->estfis = ee + ef;           // copie dans le ntuple   
      
@@ -355 +386 @@
     // C   Attention aux parametres, ici 0=OPTSHP, NO microscopic correct. 
     mglms(af,zf,0,&el);
-    // assert(isnan(el) == false);
     G4double massef = zf*fmp + (af - zf)*fmn + el + ee + ef;
-    // assert(isnan(massef) == false);
     mglms(double(aff1),double(zff1),0,&el);
-    // assert(isnan(el) == false);
     G4double masse1 = zff1*fmp + (aff1-zff1)*fmn + el + eff1;
-    // assert(isnan(masse1) == false);
     mglms(aff2,zff2,0,&el);
-    // assert(isnan(el) == false);
     G4double masse2 = zff2*fmp + (aff2 - zff2)*fmn + el + eff2;
-    // assert(isnan(masse2) == false);
     // C        WRITE(6,*) 'MASSEF,MASSE1,MASSE2',MASSEF,MASSE1,MASSE2     
     G4double b = massef - masse1 - masse2;
@@ -376 +401 @@
     } //endif
     G4double t1 = b*(b + 2.0*masse2)/(2.0*massef);
-    // assert(isnan(t1) == false);
     G4double p1 = std::sqrt(t1*(t1 + 2.0*masse1));
-    // assert(isnan(p1) == false);
     
     G4double rndm;
-    standardRandom(&rndm, &(hazard->igraine[13]));
+    rndm = randomGenerator->getRandom();
+    //    standardRandom(&rndm, &(hazard->igraine[13]));
     ctet1 = 2.0*rndm - 1.0;
-    standardRandom(&rndm,&(hazard->igraine[9]));
+    //    standardRandom(&rndm,&(hazard->igraine[9]));
+    rndm = randomGenerator->getRandom();
     phi1 = rndm*2.0*3.141592654;
            
@@ -389 +414 @@
     G4double peva = std::pow(pxeva,2) + std::pow(pyeva,2) + std::pow(pleva,2);
     G4double gamfis = std::sqrt(std::pow(massef,2) + peva)/massef;
-    // assert(isnan(gamfis) == false);
     peva = std::sqrt(peva);
-    // assert(isnan(peva) == false);
     G4double etfis = peva/massef;
       
-    G4double epf1_in;
-    G4double epf1_out;
+    G4double epf1_in = 0.0;
+    G4double epf1_out = 0.0;
 
     // C ----Matrice de rotation (noyau de fission -> Remnant)
     if(peva >= 1.0e-4) {
       sitet = std::sqrt(std::pow(pxeva,2)+std::pow(pyeva,2))/peva;
-      // assert(isnan(sitet) == false);
     }
     if(sitet > 1.0e-5) { //then
@@ -448 +470 @@
               &ffpxeva1, &ffpyeva1, &ff1, &ftype1, &inum);
       // C On ajoute le fragment:
-      // assert(af1 > 0);
       volant->iv = volant->iv + 1;
-      // assert(af1 != 0);
-      // assert(zf1 != 0);
       volant->acv[volant->iv] = af1;
       volant->zpcv[volant->iv] = zf1;
+      if(verboseLevel > 2) {
+        G4cout << __FILE__ << ":" << __LINE__ << " Added: zf1 = " << zf1 << " af1 = " << af1 << " at index " << volant->iv << G4endl;
+        volant->dump();
+      }
       if(verboseLevel > 2) {
         G4cout <<"Added fission fragment: a = " << volant->acv[volant->iv] << " z = " << volant->zpcv[volant->iv] << G4endl;
       }
       peva = std::sqrt(std::pow(ffpxeva1,2) + std::pow(ffpyeva1,2) + std::pow(ffpleva1,2));
-      // assert(isnan(peva) == false);
       volant->pcv[volant->iv] = peva;
       if(peva > 0.001) { // then
@@ -480 +502 @@
         mglms(volant->acv[iloc], volant->zpcv[iloc],0,&el);
         masse = volant->zpcv[iloc]*fmp + (volant->acv[iloc] - volant->zpcv[iloc])*fmn + el; 
-        // assert(isnan(masse) == false);
         bil1_e = bil1_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2));
-        // assert(isnan(bil1_e) == false);
         bil1_px = bil1_px + volant->pcv[iloc]*(volant->xcv[iloc]);
         bil1_py = bil1_py + volant->pcv[iloc]*(volant->ycv[iloc]);
@@ -493 +513 @@
 
       // calcul des impulsions des particules evaporees dans le systeme Remnant:
-      if(verboseLevel > 2) {
+      if(verboseLevel > 2)
         G4cout <<"2nd Translab (pf1 evap): Adding indices from " << nbpevap+1 << " to " << volant->iv << G4endl;
-      }
       nopart = varntp->ntrack - 1;
       translab(gam1,eta1,csdir1,nopart,nbpevap+1);
@@ -531 +550 @@
       volant->acv[volant->iv] = af2;
       volant->zpcv[volant->iv] = zf2; 
+      if(verboseLevel > 2)
+        G4cout << __FILE__ << ":" << __LINE__ << " Added: zf2 = " << zf2 << " af2 = " << af2 << " at index " << volant->iv << G4endl;
       if(verboseLevel > 2) {
         G4cout <<"Added fission fragment: a = " << volant->acv[volant->iv] << " z = " << volant->zpcv[volant->iv] << G4endl;
       }
       peva = std::sqrt(std::pow(ffpxeva2,2) + std::pow(ffpyeva2,2) + std::pow(ffpleva2,2));
-      // assert(isnan(peva) == false);
       volant->pcv[volant->iv] = peva;
       //      exit(0);
@@ -558 +578 @@
         masse = volant->zpcv[iloc]*fmp + (volant->acv[iloc] - volant->zpcv[iloc])*fmn + el; 
         bil2_e = bil2_e + std::sqrt(std::pow(volant->pcv[iloc],2) + std::pow(masse,2));
-        // assert(isnan(bil2_e) == false);
         bil2_px = bil2_px + volant->pcv[iloc]*(volant->xcv[iloc]);
         bil2_py = bil2_py + volant->pcv[iloc]*(volant->ycv[iloc]);
@@ -569 +588 @@
       //      G4double ctet2 = -ctet1;
       ctet2 = -1.0*ctet1;
-      assert(std::fabs(ctet2) <= 1.0);
-      // assert(isnan(ctet2) == false);
       phi2 = dmod(phi1+3.141592654,6.283185308);
-      // assert(isnan(phi2) == false);
       G4double p2 = std::sqrt(t2*(t2+2.0*masse2));
-      // assert(isnan(p2) == false);
       
       //   void translabpf(G4double masse1, G4double t1, G4double p1, G4double ctet1,
@@ -583 +598 @@
       // C  calcul des impulsions des particules evaporees dans le systeme Remnant:
       // c
-      if(verboseLevel > 2) {
-        G4cout <<"3rd Translab (pf2 evap): Adding indices from " << nbpevap+1 << " to " << volant->iv << G4endl;
-      }
+    if(verboseLevel > 2)
+      G4cout <<"3rd Translab (pf2 evap): Adding indices from " << nbpevap+1 << " to " << volant->iv << G4endl;
       nopart = varntp->ntrack - 1;
       translab(gam2,eta2,csdir2,nopart,nbpevap+1);
@@ -603 +617 @@
       
     stet1 = std::sqrt(1.0 - std::pow(ctet1,2));
-    // assert(isnan(stet1) == false);
     pf1_rem[1] = p1*stet1*std::cos(phi1);
     pf1_rem[2] = p1*stet1*std::sin(phi1);
@@ -612 +625 @@
 
     stet2 = std::sqrt(1.0 - std::pow(ctet2,2));
-    assert(std::pow(ctet2,2) >= 0.0);
-    assert(std::pow(ctet2,2) <= 1.0);
-    // assert(isnan(stet2) == false);
     
     G4double pf2_rem[4];
@@ -672 +682 @@
     // C
     //    G4double mempaw, memiv;
-    if(verboseLevel > 2) {
+    if(verboseLevel > 2)
       G4cout <<"4th Translab: Adding indices from " << memiv << " to " << volant->iv << G4endl;
-    }
     translab(gamrem,etrem,csrem,mempaw,memiv);
     // C *******************  END of fission calculations ************************
+    if(verboseLevel > 2) {
+      G4cout <<"Dump at the end of fission event " << G4endl;
+      volant->dump();
+      G4cout <<"End of dump." << G4endl;
+    }
   }
   else {
@@ -689 +703 @@
     volant->acv[volant->iv] = af;
     volant->zpcv[volant->iv] = zf;
+    if(verboseLevel > 2)
+      G4cout << __FILE__ << ":" << __LINE__ << " Added: zf = " << zf << " af = " << af << " at index " << volant->iv << G4endl;
     G4double peva = std::sqrt(std::pow(pxeva,2)+std::pow(pyeva,2)+std::pow(pleva,2));
-    // assert(isnan(peva) == false);
     volant->pcv[volant->iv] = peva;
     if(peva > 0.001) { //then
@@ -726 +741 @@
       }
       if(volant->acv[j] > 0) {
-        assert(volant->acv[j] != 0);
-        //      assert(volant->zpcv[j] != 0);
         mglms(volant->acv[j],volant->zpcv[j],0,&el);
         fmcv = volant->zpcv[j]*fmp + (volant->acv[j] - volant->zpcv[j])*fmn + el;
         e_evapo = e_evapo + std::sqrt(std::pow(volant->pcv[j],2) + std::pow(fmcv,2));
-        // assert(isnan(e_evapo) == false);
       }
     } // enddo
@@ -741 +753 @@
       
     G4double gamrem = std::sqrt(std::pow(pcorem,2)+std::pow(remmass,2))/remmass;
-    // assert(isnan(gamrem) == false);
     G4double etrem = pcorem/remmass;
 
-    if(verboseLevel > 2) {
+    if(verboseLevel > 2)
       G4cout <<"5th Translab (no fission): Adding indices from " << 1 << " to " << volant->iv << G4endl;
-    }
     nopart = varntp->ntrack - 1;
     translab(gamrem,etrem,csrem,nopart,1);
                   
+    if(verboseLevel > 2) {
+      G4cout <<"Dump at the end of evaporation event " << G4endl;
+      volant->dump();
+      G4cout <<"End of dump." << G4endl;
+    }
     // C End of the (FISSION - NO FISSION) condition (FF=1 or 0)                                          
   } //end if 
@@ -758 +773 @@
 void G4Abla::initEvapora()
 {
-  //      6     C *******************************************************************                                                                      
-  //      7     C                                                                       
-  //      8     C      SUBROUTINE ABLAINIT(STATUS,TSTAT,NAME,FPATH)
-  //      9     C********************************************************************                      
-  //     10     
-  //     11           SUBROUTINE INIT_EVAPORA(RACINE)
-  //     12                           
-  //     13     C********************************************************************                                                                       
-  //     14     C     ON INPUT:  INPUT PARAMETERS FROM FILE                             
-  //     15     C---------------------------------------------------------------------  
-  //     16     C     ON OUTPUT:                                                        
-  //     17     C     STATUS - FLAG FOR END OF INPUT FILE                               
-  //     18     C     TSTAT  - FLAG FOR NTUPLE-OUTPUT                                   
-  //     19     C     NAME   - NAME FOR ISOTOPIC PRODUCTION CROSS SECTION FILES         
-  //     20     C     FPATH  - PATH FOR  "          "        "        "    "            
-  //     21     C---------------------------------------------------------------------
-  //     22     C
-  //     23     C     Modification 5-january-2000 by KHS and BJ
-  //     24     C
-  //     25     C     New treatment of dissipation. 
-  //     26     C     See report of Beatriz Jurado, Jan. 2000
-  //     27     C
-  //     28     C---------------------------------------------------------------------
-  //     29     C   
-  //     30     C     MODIFICATION 6-aug-1999 by JB and MVR
-  //     31     C
-  //     32     C Some problems arised from an uncorrect evaluation of the fission barrier
-  //     33     C { 1) shell correction ( ECGNZ(J,K) ) was not subctracted (20-jul-99)
-  //     34     C   2) fiss. barrier EF was calc. before ang. mom. correct. (6-aug-99) }
-  //     35     C
-  //     36     C---------------------------------------------------------------------  
   //     37     C     PROJECTILE AND TARGET PARAMETERS + CROSS SECTIONS                 
   //     38     C     COMMON /ABLAMAIN/ AP,ZP,AT,ZT,EAP,BETA,BMAXNUC,CRTOT,CRNUC,       
@@ -996 +980 @@
   }
   
-  for(int z = 0; z < 98; z++) { //do 30  z = 0,98,1                                                 
+  for(int z = 0; z < 99; z++) { //do 30  z = 0,98,1                                                 
     for(int n = 0; n < 154; n++) { //do 31  n = 0,153,1                                              
       ecld->ecfnz[n][z] = 0.e0;
@@ -1003 +987 @@
       ecld->alpha[n][z] = dataInterface->getAlpha(n,z);
       ecld->vgsld[n][z] = dataInterface->getVgsld(n,z);
+      //      if(ecld->ecgnz[n][z] != 0.0) G4cout <<"ecgnz[" << n << "][" << z << "] = " << ecld->ecgnz[n][z] << G4endl;
     }
   }
@@ -1017 +1002 @@
 void G4Abla::qrot(G4double z, G4double a, G4double bet, G4double sig, G4double u, G4double *qr)
 {
-  // QROT INCLUDING DAMPING                                                
-  // INPUT: Z,A,BET,SIG,U                                                  
-  // OUTPUT: QR - COLLECTIVE ENHANCEMENT FACTOR                            
-  // 
-  // SEE  JUNGHANS ET AL., NUCL. PHYS. A 629 (1998) 635                    
-  // 
-  // 
-  // FR(U) EXPONENTIAL FUNCTION TO DEFINE DAMPING                        
-  // UCR   CRITICAL ENERGY FOR DAMPING                                   
-  // DCR   WIDTH OF DAMPING                                              
-  // BET   BETA-DEFORMATION !                                            
-  // SIG   PERPENDICULAR SPIN CUTOFF FACTOR                              
-  //   U   ENERGY                                                        
-  //  QR   COEFFICIENT OF COLLECTIVE ENHANCEMENT                         
-  //   A   MASS NUMBER                                                   
-  //   Z   CHARGE NUMBER                                                 
-
-  G4double ucr,dcr,ponq,dn,n,dz;
-
-  dcr = 10.0;
-
-  ucr = 40.0;
+  G4double ucr = 10.0; // Critical energy for damping.
+  G4double dcr = 40.0; // Width of damping.
+  G4double ponq = 0.0, dn = 0.0, n = 0.0, dz = 0.0;
 
   if(((std::fabs(bet)-1.15) < 0) || ((std::fabs(bet)-1.15) == 0)) {
@@ -1085 +1051 @@
   // USUALLY AN OBSOLETE OPTION
 
-  G4int a1,z1;
+  G4int a1 = 0, z1 = 0;
   G4double xv = 0.0, xs = 0.0, xc = 0.0, xa = 0.0;                                   
 
@@ -1143 +1109 @@
     // binding energy incl. pairing contr. is calculated from                
     // function eflmac                                                       
-    assert(a1 != 0);
     (*el) = eflmac(a1,z1,0,refopt4);
-    // assert(isnan((*el)) == false);
     if (refopt4 > 0) {
       if (refopt4 != 2) {
         (*el) = (*el) + ec2sub->ecnz[a1-z1][z1];
         //(*el) = (*el) + ec2sub->ecnz[z1][a1-z1];
-        //assert(isnan((*el)) == false);
       }
     }
@@ -1167 +1130 @@
   // RANGING FROM FISSILITY X=0.30 TO X=1.00 IN STEPS OF 0.02              
 
-  G4int index;
-  G4double x,v,dx;
+  G4int index = 0;
+  G4double x = 0.0, v = 0.0, dx = 0.0;
 
   const G4int alpha2Size = 37;
@@ -1216 +1179 @@
   //           2 : ANDREYEV                                      
 
-  G4double aa,zz,i;
+  G4double aa = 0.0, zz = 0.0, i = 0.0;
   G4double fissilityResult = 0.0;
 
@@ -1256 +1219 @@
   G4int inum = (*inum_par);
 
+  G4int idebug = 0;
+
+  if(idebug == 666) {
+    zprf =   81.;
+    aprf =   201.;
+    //    ee =   86.5877686;
+    ee = 300.0;
+    jprf =   32.;
+    zf =   0.;
+    af =   0.;
+    mtota =   0.;
+    pleva =   0.;
+    pxeva =   0.;
+    pyeva =   0.;
+    ff =  -1;
+    inttype =  0;
+    inum =  2;
+    G4cout <<" PK::: EVAPORA event  " << inum << G4endl;
+    G4cout <<" PK:::    zprf =   " << zprf << G4endl;
+    G4cout <<" PK:::    aprf =   " << aprf << G4endl;
+    G4cout <<" PK:::      ee =   " << ee << G4endl;
+    G4cout <<" PK:::      eeDiff =   " << ee - 86.5877686 << G4endl;
+    G4cout <<" PK:::    jprf =   " << jprf << G4endl;
+    G4cout <<" PK:::      zf =   " << zf << G4endl;
+    G4cout <<" PK:::      af =   " << af << G4endl;
+    G4cout <<" PK:::   mtota =   " << mtota << G4endl;
+    G4cout <<" PK:::   pleva =   " << pleva << G4endl;
+    G4cout <<" PK:::   pxeva =   " << pxeva << G4endl;
+    G4cout <<" PK:::   pyeva =   " << pyeva << G4endl;
+    G4cout <<" PK:::      ff =  " << ff << G4endl;
+    G4cout <<" PK::: inttype =  " << inttype << G4endl;
+    G4cout <<" PK:::    inum =  " << inum << G4endl;
+  }
   //    533     C                                                                       
   //    534     C     INPUT:                                                            
@@ -1348 +1344 @@
   // SAVE -> static
         
-  static G4int sortie;                            
-  static G4double epsiln,probp,probn,proba,ptotl,e;  
-  static G4double sn,sbp,sba,x,amoins,zmoins,ecn,ecp,eca,bp,ba;         
-  static G4double pteva;                       
-
-  static G4int itest;
-  static G4double probf;
-
-  static G4int k, j, il;
-
-  static G4double ctet1,stet1,phi1;
-  static G4double sbfis,rnd;
-  static G4double selmax;
-  static G4double segs;
-  static G4double ef;
-  static G4int irndm;
-
-  static G4double pc, malpha;
+  static G4int sortie = 0;                            
+  static G4double epsiln = 0.0, probp = 0.0, probn = 0.0, proba = 0.0, ptotl = 0.0, e = 0.0;  
+  static G4double sn = 0.0, sbp = 0.0, sba = 0.0, x = 0.0, amoins = 0.0, zmoins = 0.0;
+  G4double ecn = 0.0, ecp = 0.0,eca = 0.0, bp = 0.0, ba = 0.0;         
+  static G4double pteva = 0.0;                       
+
+  static G4int itest = 0;
+  static G4double probf = 0.0;
+
+  static G4int k = 0, j = 0, il = 0;
+
+  static G4double ctet1 = 0.0, stet1 = 0.0, phi1 = 0.0;
+  static G4double sbfis = 0.0, rnd = 0.0;
+  static G4double selmax = 0.0;
+  static G4double segs = 0.0;
+  static G4double ef = 0.0;
+  static G4int irndm = 0;
+
+  static G4double pc = 0.0, malpha = 0.0;
 
   zf = zprf;
@@ -1392 +1389 @@
   direct(zf,af,ee,jprf,&probp,&probn,&proba,&probf,&ptotl,
          &sn,&sbp,&sba,&ecn,&ecp,&eca,&bp,&ba,inttype,inum,itest); //:::FIXME::: Call
-  // assert(isnan(proba) == false);
-  // assert(isnan(probp) == false);  
-  // assert(isnan(probn) == false);
-  // assert(isnan(probf) == false);  
-  assert((eca+ba) >= 0);
-  assert((ecp+bp) >= 0);
-  // assert(isnan(ecp) == false);
-  // assert(isnan(ecn) == false);
-  // assert(isnan(bp) == false);
-  // assert(isnan(ba) == false);
+  if((eca+ba) < 0) {
+    eca = 0.0;
+    ba = 0.0;
+  }
   k = idnint(zf);
   j = idnint(af-zf);
@@ -1411 +1402 @@
   il  = idnint(jprf);
   barfit(k,k+j,il,&sbfis,&segs,&selmax);
-  // assert(isnan(sbfis) == false);
   
   if ((fiss->optshp == 1) || (fiss->optshp == 3)) { //then                     
@@ -1423 +1413 @@
   ef = fb->efa[j][k];
   //  ef = fb->efa[j][k];
-  // assert(isnan(fb->efa[j][k]) == false);
   // here the final steps of the evaporation are calculated                
   if ((sortie == 1) || (ptotl == 0.e0)) {
@@ -1460 +1449 @@
   // random number for the evaporation                                     
   //  x = double(Rndm(irndm))*ptotl;
-  x = double(haz(1))*ptotl;
-
+  // x = double(haz(1))*ptotl;
+  x = randomGenerator->getRandom() * ptotl;
 //   G4cout <<"proba = " << proba << G4endl;
 //   G4cout <<"probp = " << probp << G4endl;
@@ -1471 +1460 @@
     // alpha evaporation                                                     
     if (itest == 1) {
-      G4cout <<"< alpha evaporation >" << G4endl;
+      G4cout <<"PK::: < alpha evaporation >" << G4endl;
     }
     amoins = 4.0;
     zmoins = 2.0;
     epsiln = sba + eca;
-    assert((std::pow((1.0 + (eca+ba)/3.72834e3),2) - 1.0) >= 0);
     pc = std::sqrt(std::pow((1.0 + (eca+ba)/3.72834e3),2) - 1.0) * 3.72834e3;
-    // assert(isnan(pc) == false);
     malpha = 4.0;
 
@@ -1490 +1477 @@
     // proton evaporation                                                    
     if (itest == 1) {
-      G4cout <<"< proton evaporation >" << G4endl;
+      G4cout <<"PK::: < proton evaporation >" << G4endl;
     }
     amoins = 1.0;
     zmoins = 1.0;
     epsiln = sbp + ecp;
-    assert((std::pow((1.0 + (ecp + bp)/9.3827e2),2) - 1.0) >= 0);
     pc = std::sqrt(std::pow((1.0 + (ecp + bp)/9.3827e2),2) - 1.0) * 9.3827e2;
-    // assert(isnan(pc) == false);
     malpha = 0.0;
     // volant:
     volant->iv = volant->iv + 1;
-    volant->acv[volant->iv] = 1.;
+    volant->acv[volant->iv] = 1.0;
     volant->zpcv[volant->iv] = 1.;
     volant->pcv[volant->iv] = pc;
@@ -1508 +1493 @@
     // neutron evaporation                                                   
     if (itest == 1) {
-      G4cout <<"< neutron evaporation >" << G4endl;
+      G4cout <<"PK::: < neutron evaporation >" << G4endl;
     }
     amoins = 1.0;
     zmoins = 0.0;
     epsiln = sn + ecn;
-    assert((std::pow((1.0 + (ecn)/9.3956e2),2) - 1.0) >= 0);
     pc = std::sqrt(std::pow((1.0 + (ecn)/9.3956e2),2) - 1.0) * 9.3956e2;
-    // assert(isnan(pc) == false);
+    if(itest == 1) {
+      G4cout <<"PK::: pc " << pc << G4endl;
+    }
     malpha = 0.0;
   
@@ -1523 +1509 @@
     volant->zpcv[volant->iv] = 0.;
     volant->pcv[volant->iv] = pc;
+
+    if(volant->getTotalMass() > 209 && verboseLevel > 0) {
+      volant->dump();
+      G4cout <<"DEBUGA Total = " << volant->getTotalMass() << G4endl;
+    }
   }
   else {
@@ -1530 +1521 @@
     // fission fragment mass distribution is calulated in subroutine fisdis  
     if (itest == 1) {
-      G4cout <<"< fission >" << G4endl;
+      G4cout <<"PK::: < fission >" << G4endl;
     }
     amoins = 0.0;
@@ -1539 +1530 @@
     pc = 0.0;
     ff = 1;
-    //    ff = 0; // For testing, allows to disable fission!
+    // ff = 0; // For testing, allows to disable fission!
   }
 
   if (itest == 1) {
-    G4cout <<"sn,sbp,sba,ef" << sn << "," << sbp << "," << sba <<"," << ef << G4endl;
-    G4cout <<"probn,probp,proba,probf,ptotl " <<","<< probn <<","<< probp <<","<< proba <<","<< probf <<","<< ptotl << G4endl;
+    G4cout << std::setprecision(9) <<"PK::: SN,SBP,SBA,EF  " << sn << "  " << sbp << "  " << sba <<"  " << ef << G4endl;
+    G4cout << std::setprecision(9) <<"PK::: PROBN,PROBP,PROBA,PROBF,PTOTL  " <<"  "<< probn <<"  "<< probp <<"  "<< proba <<"  "<< probf <<"  "<< ptotl << G4endl;
   }
 
@@ -1557 +1548 @@
 
   if(ff == 0) {
-    standardRandom(&rnd,&(hazard->igraine[8]));
+    rnd = randomGenerator->getRandom();
+    //    standardRandom(&rnd,&(hazard->igraine[8]));
     ctet1 = 2.0*rnd - 1.0;
-    standardRandom(&rnd,&(hazard->igraine[4]));
+    rnd = randomGenerator->getRandom();
+    //    standardRandom(&rnd,&(hazard->igraine[4]));
     phi1 = rnd*2.0*3.141592654;
     stet1 = std::sqrt(1.0 - std::pow(ctet1,2));
-    // assert(isnan(stet1) == false);
     volant->xcv[volant->iv] = stet1*std::cos(phi1);
     volant->ycv[volant->iv] = stet1*std::sin(phi1);
@@ -1569 +1561 @@
     pyeva = pyeva - pc * volant->ycv[volant->iv];
     pleva = pleva - pc * ctet1;
-    // assert(isnan(pleva) == false);
   }
 
@@ -1599 +1590 @@
                     G4double *ba_par, G4int inttype, G4int inum, G4int itest)
 {
-  G4int dummy0;
+  G4int dummy0 = 0;
   
   G4double probp = (*probp_par);
@@ -1680 +1671 @@
   // - AVERAGE KINETIC ENERGIES (2*T) AND EFFECTIVE BARRIERS     
 
-  static G4double bk;
-  static G4int afp;
-  static G4double at; // = 0.0;
-  static G4double bet; // = 0.0;
-  static G4double bs;
-  static G4double bshell;
-  static G4double cf;
-  static G4double dconst;
-  static G4double defbet;
-  static G4double denomi;
-  static G4double densa;
-  static G4double densf;
-  static G4double densg;
-  static G4double densn;
-  static G4double densp;
-  static G4double edyn;
-  static G4double eer;
-  static G4double ef;
-  static G4double ft;
-  static G4double ga; // = 0.0;
-  static G4double gf; // = 0.0;
-  static G4double gn; // = 0.0;
-  static G4double gngf;
-  static G4double gp; // = 0.0;
-  static G4double gsum;
+  static G4double bk = 0.0;
+  static G4int afp = 0;
+  static G4double at = 0.0;
+  static G4double bs = 0.0;
+  static G4double bshell = 0.0;
+  static G4double cf = 0.0;
+  static G4double dconst = 0.0;
+  static G4double defbet = 0.0;
+  static G4double denomi = 0.0;
+  static G4double densa = 0.0;
+  static G4double densf = 0.0;
+  static G4double densg = 0.0;
+  static G4double densn = 0.0;
+  static G4double densp = 0.0;
+  static G4double edyn = 0.0;
+  static G4double eer = 0.0;
+  static G4double ef = 0.0;
+  static G4double ft = 0.0;
+  static G4double ga = 0.0;
+  static G4double gf = 0.0;
+  static G4double gn = 0.0;
+  static G4double gngf = 0.0;
+  static G4double gp = 0.0;
+  static G4double gsum = 0.0;
   static G4double hbar = 6.582122e-22; // = 0.0;
-  static G4double homega; // = 0.0;
-  static G4double iflag;
-  static G4int il;
-  static G4int ilast;
-  static G4int imaxwell;
-  static G4int in;
-  static G4int iz;
-  static G4int j;
-  static G4int k;
-  static G4double ma1z;
-  static G4double ma1z1;
-  static G4double ma4z2;
-  static G4double maz;
-  static G4double nprf;
-  static G4double nt;
-  static G4double parc;
+  static G4double iflag = 0.0;
+  static G4int il = 0;
+  static G4int imaxwell = 0;
+  static G4int in = 0;
+  static G4int iz = 0;
+  static G4int j = 0;
+  static G4int k = 0;
+  static G4double ma1z = 0.0;
+  static G4double ma1z1 = 0.0;
+  static G4double ma4z2 = 0.0;
+  static G4double maz = 0.0;
+  static G4double nprf = 0.0;
+  static G4double nt = 0.0;
+  static G4double parc = 0.0;
   static G4double pi = 3.14159265;
-  static G4double pt;
-  static G4double ra;
-  static G4double rat;
-  static G4double refmod;
-  static G4double rf;
-  static G4double rn;
-  static G4double rnd;
-  static G4double rnt;
-  static G4double rp;
-  static G4double rpt;
-  static G4double sa;
-  static G4double sbf;
-  static G4double sbfis;
-  static G4double segs;
-  static G4double selmax;
-  static G4double sp;
-  static G4double tauc;
-  static G4double tconst;
-  static G4double temp;
-  static G4double ts1;
-  static G4double tsum;
-  static G4double wf; // = 0.0;
-  static G4double wfex;
-  static G4double xx;
-  static G4double y;
+  static G4double pt = 0.0;
+  static G4double ra = 0.0;
+  static G4double rat = 0.0;
+  static G4double refmod = 0.0;
+  static G4double rf = 0.0;
+  static G4double rn = 0.0;
+  static G4double rnd = 0.0;
+  static G4double rnt = 0.0;
+  static G4double rp = 0.0;
+  static G4double rpt = 0.0;
+  static G4double sa = 0.0;
+  static G4double sbf = 0.0;
+  static G4double sbfis = 0.0;
+  static G4double segs = 0.0;
+  static G4double selmax = 0.0;
+  static G4double sp = 0.0;
+  static G4double tauc = 0.0;
+  static G4double tconst = 0.0;
+  static G4double temp = 0.0;
+  static G4double ts1 = 0.0;
+  static G4double tsum = 0.0;
+  static G4double wf = 0.0;
+  static G4double wfex = 0.0;
+  static G4double xx = 0.0;
+  static G4double y = 0.0;
 
   imaxwell = 1;
@@ -1757 +1745 @@
 
   // no limit if statistical model is calculated.                         
-  if (bet <= 1.0e-16) {
+  if (fiss->bet <= 1.0e-16) {
     edyn = 10000.0;
   }
@@ -1763 +1751 @@
   // just a change of name until the end of this subroutine                
   eer = ee;
+  if(verboseLevel > 2)
+    G4cout << __FILE__ << ":" << __LINE__ << " eer = " << eer << G4endl;
   if (inum == 1) {
     ilast = 1;
@@ -1784 +1774 @@
     mglw(a-4.0,zprf-2.0,&ma4z2);
   }
-  // assert(isnan(maz) == false);
-  // assert(isnan(ma1z) == false);
-  // assert(isnan(ma1z1) == false);
-  // assert(isnan(ma4z2) == false);
   
   // separation energies and effective barriers                     
@@ -1804 +1790 @@
 
   sbp = sp + bp;
-  // assert(isnan(sbp) == false);
-  // assert(isinf(sbp) == false);
   if (a-4.0 <= 0.0) {
     sba = 1.0e+75;
@@ -1817 +1801 @@
 
   sba = sa + ba;
-  // assert(isnan(sba) == false);
-  // assert(isinf(sba) == false);
  direct30:
 
@@ -1833 +1815 @@
     il = idnint(jprf);
     barfit(k,k+j,il,&sbfis,&segs,&selmax);
-    // assert(isnan(sbfis) == false);
     if ((fiss->optshp == 1) || (fiss->optshp == 3)) {
       //      fb->efa[k][j] = G4double(sbfis) -  ecld->ecgnz[j][k];
@@ -1857 +1838 @@
       fb->efa[j][k] = 0.0;
     }
-    // assert(isnan(fb->efa[j][k]) == false);
     
     // factor with jprf should be 0.0025d0 - 0.01d0 for                     
@@ -1867 +1847 @@
     }
     xx = fissility((k+j),k,fiss->optxfis);
-    // assert(isnan(xx) == false);
-    // assert(isinf(xx) == false);
     
     y = 1.00 - xx;
@@ -1878 +1856 @@
     }
     bs = bipol(1,y);
-    // assert(isnan(bs) == false);
-    // assert(isinf(bs) == false);
     bk = bipol(2,y);
-    // assert(isnan(bk) == false);
-    // assert(isinf(bk) == false);
   }
   else {
@@ -1895 +1869 @@
   in = afp - iz;
   bshell = ecld->ecfnz[in][iz];
-  // assert(isnan(bshell) == false);
 
   // ld saddle point deformation                                          
@@ -1905 +1878 @@
   // a.r.j. 6.8.97                                                         
   defbet = 1.58533e0 * spdef(idnint(a),idnint(zprf),fiss->optxfis);
-  // assert(isnan(defbet) == false);
   
   // level density and temperature at the saddle point                     
@@ -1919 +1891 @@
   //   G4cout <<"densf = " << densf << G4endl;
   //   G4cout <<"temp = " << temp << G4endl;
-  // assert(isnan(densf) == false);
-  // assert(isnan(temp) == false);
-  //  assert(temp != 0);
   ft = temp;
   if (iz >= 2) {
@@ -1929 +1898 @@
     // level density and temperature in the proton daughter                  
     densniv(a-1.0,zprf-1.0e0,ee,sbp,&densp, bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
-    assert(temp >= 0);
-    // assert(isnan(temp) == false);
     pt = temp;
     if (imaxwell == 1) {
@@ -1946 +1913 @@
       iflag = iflag + 1;
       if(iflag >= 10) {
-        standardRandom(&rnd,&(hazard->igraine[5]));
+        rnd = randomGenerator->getRandom();
+        //      standardRandom(&rnd,&(hazard->igraine[5]));
         ecp=std::sqrt(rnd)*(eer-sbp);
-        // assert(isnan(ecp) == false);
         goto direct2914;
       }
@@ -1976 +1943 @@
     densniv(a-1.0,zprf,ee,sn,&densn,bshell, 1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
     nt = temp;
-
+    if(itest == 1)
+      G4cout <<"PK::: nt = " << nt <<G4endl;
     if (imaxwell == 1) {
       // valentina - random kinetic energy in a maxwelliam distribution
@@ -1988 +1956 @@
 
       iflag=0;
-
+    direct1915:
       ecn = fmaxhaz(rnt);
+      if(verboseLevel > 2) {
+        G4cout <<"rnt = " << rnt << G4endl;
+        G4cout << __FILE__ << ":" << __LINE__ << " ecn = " << ecn << G4endl;
+      }
       iflag=iflag+1;
       if(iflag >= 10) {
-        standardRandom(&rnd,&(hazard->igraine[6]));
+        rnd = randomGenerator->getRandom();
+        //      standardRandom(&rnd,&(hazard->igraine[6]));
         ecn = std::sqrt(rnd)*(eer-sn);
-        // assert(isnan(ecn) == false);
+        if(verboseLevel > 2)
+          G4cout << __FILE__ << ":" << __LINE__ << " ecn = " << ecn << G4endl;
         goto direct2915;
       }
-      //       if((ecn+sn) > eer) {
-      //        goto direct1915;
-      //       }
-      //       else {
-      //        ecn = 2.e0 * nt;
-      //       }
-      if((ecn + sn) <= eer) {
-        ecn = 2.0 * nt;
-      }
+      if((ecn+sn) > eer) {
+        goto direct1915;
+      } 
+    } 
+    else {
+        ecn = 2.e0 * nt;
+        if(verboseLevel > 2)
+          G4cout << __FILE__ << ":" << __LINE__ << " ecn = " << ecn << G4endl;
+    }
+//       if((ecn + sn) <= eer) {
+//      ecn = 2.0 * nt;
+//      G4cout << __FILE__ << ":" << __LINE__ << " ecn = " << ecn << G4endl;
+//       }
     direct2915: 
       dummy0 = 0;
       //      G4cout <<"" <<G4endl;
-    }
-  }
+  } 
   else {
     densn = 0.0;
@@ -2021 +1998 @@
     defbet = 1.5 * (ecld->alpha[in-2][iz-2]);
 
+    // For debugging:
+    //    bshell = -10.7;
+    //    defbet = -0.06105;
+    // G4cout <<"ecgnz N = " << in-2 << G4endl;
+    // G4cout <<"ecgnz Z = " << iz-2 << G4endl;
+    //  G4cout <<"bshell = " << bshell << G4endl;
+    // G4cout <<"defbet = " << defbet << G4endl;
     // level density and temperature in the alpha daughter                   
     densniv(a-4.0,zprf-2.0e0,ee,sba,&densa,bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
+    // G4cout <<"densa = " << densa << G4endl;
+    // G4cout <<"temp = " << temp << G4endl;
 
     // valentina - random kinetic energy in a maxwelliam distribution
@@ -2039 +2025 @@
       iflag=iflag+1;
       if(iflag >= 10) {
-        standardRandom(&rnd,&(hazard->igraine[7]));
+        rnd = randomGenerator->getRandom();
+        //      standardRandom(&rnd,&(hazard->igraine[7]));
         eca=std::sqrt(rnd)*(eer-sba);
-        // assert(isnan(eca) == false);
         goto direct2916;
       }
@@ -2047 +2033 @@
         goto direct1916;
       }
-      else {
-        eca = 2.0 * at;
-      }
+    }
+    else {
+      eca = 2.0 * at;
+    }
     direct2916:
       dummy0 = 0;
       //      G4cout <<"" << G4endl;
-    }
-    else {
-      densa = 0.0;
-      eca = 0.0;
-      at = 0.0;
-    }
-  } // PK
+  }
+  else {
+    densa = 0.0;
+    eca = 0.0;
+    at = 0.0;
+  }
+  //} // PK
 
   // special treatment for unbound nuclei                                                
@@ -2087 +2074 @@
   // compound nucleus level density                                        
   densniv(a,zprf,ee,0.0e0,&densg,bshell,1.e0,1.e0,&temp,int(fiss->optshp),int(fiss->optcol),defbet);
-  // assert(isnan(densg) == false);
-  // assert(isnan(temp) == false);
   
   if ( densg > 0.e0) {
@@ -2097 +2082 @@
     ga = (std::pow(a,(2.0/3.0))/fiss->akap)*densa/densg/pi*2.0*std::pow(at,2);
     gf = densf/densg/pi/2.0*ft;
-    // assert(isnan(gf) == false);
-    
-    //     assert(isnan(gp) == false);
-    //     assert(isnan(gn) == false);
-    //     assert(isnan(ga) == false);
-    //     assert(isnan(ft) == false);
-    //    assert(ft != 0);
-    //    assert(isnan(gf) == false);
     
     if(itest == 1) {
@@ -2117 +2094 @@
 
   gsum = ga + gp + gn;
-  // assert(isinf(gsum) == false);
-  // assert(isnan(gsum) == false);
   if (gsum > 0.0) {
     ts1  = hbar / gsum;
@@ -2162 +2137 @@
     rp = 1.0;
     rn = densn/densp*std::pow((nt/pt),2);
-    // assert(isnan(rn) == false);
     ra = densa*2.0/densp*std::pow((at/pt),2);
-    // assert(isnan(ra) == false);
     goto direct50;
   }
@@ -2173 +2146 @@
   // cramers and weidenmueller factors for the dynamical hindrances of     
   // fission                                                               
-  if (bet <= 1.0e-16) {
+  if (fiss->bet <= 1.0e-16) {
     cf = 1.0;
     wf = 1.0;
   }
   else if (sbf > 0.0e0) {
-    cf = cram(bet,homega);
+    cf = cram(fiss->bet,fiss->homega);
     // if fission barrier ef=0.d0 then fission is the only possible      
     // channel. to avoid std::log(0) in function tau                          
@@ -2190 +2163 @@
     else {
       // transient time tau()                                                  
-      tauc = tau(bet,homega,ef,ft);
-      // assert(isnan(tauc) == false);
+      tauc = tau(fiss->bet,fiss->homega,ef,ft);
     }
     wfex = (tauc - tsum)/ts1;
@@ -2216 +2188 @@
   tconst = 0.7;
   dconst = 12.0/std::sqrt(a);
-  // assert(isnan(dconst) == false);
   nprf = a - zprf;
 
   if (fiss->optshp >= 2) { //then                                           
     parite(nprf,&parc);
-    // assert(isnan(parc) == false);
     dconst = dconst*parc;
   }
@@ -2241 +2211 @@
     else {
       rn=gngf;
-      // assert(isnan(rn) == false);
       rp=gngf*gp/gn;
-      // assert(isnan(rp) == false);      
       ra=gngf*ga/gn;
-      // assert(isnan(ra) == false);      
     }
   } else {
-    // assert(isnan(cf) == false);
-    // assert(isinf(gn) == false);
-    // assert(isinf(gf) == false);
-    // assert(isinf(cf) == false);        
-    assert(gn > 0 || (gf != 0 && cf != 0));
     rn = gn/(gf*cf);
 //     G4cout <<"rn = " << G4endl;
 //     G4cout <<"gn = " << gn << " gf = " << gf << " cf = " << cf << G4endl;
-    // assert(isnan(rn) == false);
     rp = gp/(gf*cf);
-    // assert(isnan(rp) == false);    
     ra = ga/(gf*cf);
-    // assert(isnan(ra) == false);    
   }
  direct50:
   // relative decay probabilities                                          
-  // assert(isnan(ra) == false);
-  // assert(isnan(rp) == false);
-  // assert(isnan(rn) == false);
-  // assert(isnan(rf) == false);
-  
   denomi = rp+rn+ra+rf;
-  // assert(isnan(denomi) == false);
-  assert(denomi > 0);
   // decay probabilities after transient time
   probf = rf/denomi;
-  // assert(isnan(probf) == false);
   probp = rp/denomi;
-  // assert(isnan(probp) == false);
   probn = rn/denomi;
-  // assert(isnan(probn) == false);
   proba = ra/denomi;
-  // assert(isnan(proba) == false);
-  // assert(isinf(proba) == false);
   
   // new treatment of grange-weidenmueller factor, 5.1.2000, khs !!!
 
   // decay probabilites with transient time included
-  // assert(isnan(wf) == false);
-  assert(std::fabs(probf) <= 1.0);
   probf = probf * wf;
   if(probf == 1.0) {
@@ -2300 +2245 @@
   }
  direct70:
-  // assert(isnan(probp) == false);
-  // assert(isnan(probn) == false);
-  // assert(isnan(probf) == false);  
-  // assert(isnan(proba) == false);
   ptotl = probp+probn+proba+probf;
-  // assert(isnan(ptotl) == false);
   
   ee = eer;
@@ -2311 +2251 @@
 
   // Return values:
-  // assert(isnan(proba) == false);
   (*probp_par) = probp;
   (*probn_par) = probn;
@@ -2378 +2317 @@
   //   1544     C                                                                       
   //   1545     C                                                                       
-  G4double afp;
-  G4double delta0;
-  G4double deltau;
-  G4double deltpp;
-  G4double e;
+  G4double afp = 0.0;
+  G4double delta0 = 0.0;
+  G4double deltau = 0.0;
+  G4double deltpp = 0.0;
+  G4double e = 0.0;
   G4double ecor = 0.0;
-  G4double ecor1;
-  G4double ecr;
-  G4double er;
-  G4double fe;
-  G4double fp;
-  G4double he;
-  G4double iz;
-  G4double pa;
-  G4double para;
-  G4double parz;
-  G4double ponfe;
-  G4double ponniv;
-  G4double qr;
-  G4double sig;
-  G4double y01;
-  G4double y11;
-  G4double y2;
-  G4double y21;
-  G4double y1;
-  G4double y0;
+  G4double ecor1 = 0.0;
+  G4double ecr = 0.0;
+  G4double er = 0.0;
+  G4double fe = 0.0;
+  G4double fp = 0.0;
+  G4double he = 0.0;
+  G4double iz = 0.0;
+  G4double pa = 0.0;
+  G4double para = 0.0;
+  G4double parz = 0.0;
+  G4double ponfe = 0.0;
+  G4double ponniv = 0.0;
+  G4double qr = 0.0;
+  G4double sig = 0.0;
+  G4double y01 = 0.0;
+  G4double y11 = 0.0;
+  G4double y2 = 0.0;
+  G4double y21 = 0.0;
+  G4double y1 = 0.0;
+  G4double y0 = 0.0;
 
   G4double pi6 = std::pow(3.1415926535,2) / 6.0;
@@ -2453 +2392 @@
       //      deltpp = -0.25e0* (delta0/std::pow(std::sqrt(a),2)) * pa /pi6 + 2.e0*delta0/std::sqrt(a);
       deltpp = -0.25e0* std::pow((delta0/std::sqrt(a)),2) * pa /pi6 + 2.e0*delta0/std::sqrt(a);
-      // assert(isnan(deltpp) == false);
       
       parite(a,&para);
       if (para < 0.0) {
         e = e - delta0/std::sqrt(a);
-        // assert(isnan(e) == false);
       } else {                                                         
         parite(z,&parz);
         if (parz > 0.e0) {
           e = e - 2.0*delta0/std::sqrt(a);
-          // assert(isnan(e) == false);
         } else {
           e = e;
-          // assert(isnan(e) == false);
         }
       }
@@ -2510 +2445 @@
   if (ee < 5.0) {
     y1 = std::sqrt(pa*ecor);
-    // assert(isnan(y1) == false);
     for(int j = 0; j < 5; j++) {
       y2 = pa*ecor*(1.e0-std::exp(-y1));
-      // assert(isnan(y2) == false);
       y1 = std::sqrt(y2);
-      // assert(isnan(y1) == false);
     }
     
     y0 = pa/y1;
-    // assert(isnan(y0) == false);
-    assert(y0 != 0.0);
     (*temp)=1.0/y0;
     (*dens) = std::exp(y0*ecor)/ (std::pow((std::pow(ecor,3)*y0),0.5)*std::pow((1.0-0.5*y0*ecor*std::exp(-y1)),0.5))* std::exp(y1)*(1.0-std::exp(-y1))*0.1477045;
@@ -2526 +2456 @@
       ecor1=1.0;
       y11 = std::sqrt(pa*ecor1);
-      // assert(isnan(y11) == false);
       for(int j = 0; j < 7; j++) {
         y21 = pa*ecor1*(1.0-std::exp(-y11));
-        // assert(isnan(21) == false);
         y11 = std::sqrt(y21);
-        // assert(isnan(y11) == false);
       }
 
       y01 = pa/y11;
-      // assert(isnan(y01) == false);
       (*dens) = (*dens)*std::pow((y01/y0),1.5);
       (*temp) = (*temp)*std::pow((y01/y0),1.5);
@@ -2542 +2468 @@
   else {
     ponniv = 2.0*std::sqrt(pa*ecor);
-    // assert(isnan(ponniv) == false);
     if (ponniv > 700.0) {
       ponniv = 700.0;
@@ -2549 +2474 @@
     // fermi gas state density                                               
     (*dens) = std::pow(pa,(-0.25e0))*std::pow(ecor,(-1.25e0))*std::exp(ponniv) * 0.1477045e0;
-    // assert(isnan(std::sqrt(ecor/pa)) == false);
     (*temp) = std::sqrt(ecor/pa);
   }
@@ -2566 +2490 @@
 
   (*dens) = (*dens) * qr;
+  if(verboseLevel > 2) {
+    G4cout <<"PK::: dens = " << (*dens) << G4endl;
+    G4cout <<"PK::: AFP, IZ, ECOR, ECOR1 " << afp << " " << iz << " " << ecor << " " << ecor1 << G4endl;
+  }
 }
 
@@ -2576 +2504 @@
   // replaces tables. Barrier heights from Myers and Swiatecki !!!                                                                 
 
-  G4double nms,ims,ksims,xms, ums;
+  G4double nms = 0.0, ims = 0.0, ksims = 0.0, xms = 0.0, ums = 0.0;
 
   nms = ams - zms;
@@ -2623 +2551 @@
   //       otherwise   - return pairing (= -1/2 dpn + 1/2 (Dp + Dn))   
 
-  G4double eflmacResult;
-
-  G4int in;
-  G4double z,n,a,av,as,a0,c1,c4,b1,b3,f,ca,w,dp,dn,dpn,efl,pi;
-  G4double rmac,bs,h,r0,kf,ks,kv,rp,ay,aden,x0,y0,mh,mn,esq,ael,i;
-  pi = 3.141592653589793238e0;
+  G4double eflmacResult = 0.0;
+
+  G4int in = 0;
+  G4double z = 0.0, n = 0.0, a = 0.0, av = 0.0, as = 0.0;
+  G4double a0 = 0.0, c1 = 0.0, c4 = 0.0, b1 = 0.0, b3 = 0.0;
+  G4double f = 0.0, ca = 0.0, w = 0.0, dp = 0.0, dn = 0.0, dpn = 0.0, efl = 0.0;
+  G4double rmac = 0.0, bs = 0.0, h = 0.0, r0 = 0.0, kf = 0.0, ks = 0.0;
+  G4double kv = 0.0, rp = 0.0, ay = 0.0, aden = 0.0, x0 = 0.0, y0 = 0.0;
+  G4double mh = 0.0, mn = 0.0, esq = 0.0, ael = 0.0, i = 0.0;
+  G4double pi = 3.141592653589793238e0;
 
   // fundamental constants
@@ -2696 +2628 @@
   dp  = rmac*bs/std::pow(z,(1.0/3.0));
   dpn = h/bs/std::pow(a,(2.0/3.0));
-  // assert(isnan(dpn) == false);
   
   c1  = 3.0/5.0*esq/r0;
-  // assert(isnan(c1) == false);
-  // assert(isinf(c1) == false);
-  
   c4  = 5.0/4.0*std::pow((3.0/(2.0*pi)),(2.0/3.0)) * c1;
-  // assert(isnan(c4) == false);
-  // assert(isinf(c4) == false);
- 
-  // assert(isnan(pi) == false);
-  // assert(isnan(z) == false);
-  // assert(isnan(a) == false);
-  // assert(isnan(r0) == false);
   kf  = std::pow((9.0*pi*z/(4.0*a)),(1.0/3.0))/r0;
-  // assert(isnan(kf) == false);
-  // assert(isinf(kf) == false);
   
   f = -1.0/8.0*rp*rp*esq/std::pow(r0,3) * (145.0/48.0 - 327.0/2880.0*std::pow(kf,2) * std::pow(rp,2) + 1527.0/1209600.0*std::pow(kf,4) * std::pow(rp,4));
@@ -2772 +2691 @@
   // ENERGY OF A SPECIFIC NUCLEUS                                      
 
-  double para,parz;
+  double para = 0.0, parz = 0.0;
   // A                 MASS NUMBER                                     
   // Z                 NUCLEAR CHARGE                                  
@@ -2787 +2706 @@
     parite(z, &parz);
     if (parz > 0.0) {
-      // assert(isnan(std::sqrt(a)) == false);
       (*del) = -12.0/std::sqrt(a);
     }
     else {
-      // assert(isnan(std::sqrt(a)) == false);
       (*del) = 12.0/std::sqrt(a);
     }
@@ -2804 +2721 @@
   // RETURNS -1 IF N IS ODD AND +1 IF N IS EVEN                        
 
-  G4double n1, n2, n3;
+  G4double n1 = 0.0, n2 = 0.0, n3 = 0.0;
 
   // N                 NUMBER TO BE TESTED                             
@@ -2833 +2750 @@
   // T      - NUCLEAR TEMPERATURE                                        
 
-  G4double tauResult;
-
-  G4double tlim;
-  tlim = 8.e0 * ef;
+  G4double tauResult = 0.0;
+
+  G4double tlim = 8.e0 * ef;
   if (t > tlim) {
     t = tlim;
@@ -2844 +2760 @@
   if (bet/(std::sqrt(2.0)*10.0*(homega/6.582122)) <= 1.0) {
     tauResult = std::log(10.0*ef/t)/(bet*1.0e21);
-    //    assert(isnan(tauResult) == false);
   }
   else {
     tauResult = std::log(10.0*ef/t)/ (2.0*std::pow((10.0*homega/6.582122),2))*(bet*1.0e-21);
-    // assert(isnan(tauResult) == false);
   } //end if                                                            
 
@@ -2868 +2782 @@
   }
 
-  // assert(isnan(cramResult) == false);
   return cramResult;
 }
@@ -2883 +2796 @@
   // LINEAR INTERPOLATION OF BS BK TABLE                               
 
-  int i;
-
-  G4double bipolResult;
+  int i = 0;
+
+  G4double bipolResult = 0.0;
 
   const int bsbkSize = 54;
@@ -2908 +2821 @@
 
   i = idint(y/(2.0e-02)) + 1;
-  assert(i >= 1);
     
   if(i >= bsbkSize) {
@@ -2925 +2837 @@
   }
   
-  // assert(isnan(bipolResult) == false);
   return bipolResult;
 }
@@ -2998 +2909 @@
 
   G4double pa[7],pz[7],pl[10];
-  G4double a,z,amin,amax,amin2,amax2,aa,zz,bfis;
-  G4double bfis0,ell,el,egs,el80,el20,elmax,sel80,sel20,x,y,q,qa,qb;
-  G4double aj,ak,a1,a2;
-
-  G4int i,j,k,m;
-  G4int l;
+  for(G4int init_i = 0; init_i < 7; init_i++) {
+    pa[init_i] = 0.0; 
+    pz[init_i] = 0.0; 
+  }
+  for(G4int init_i = 0; init_i < 10; init_i++) {
+    pl[init_i] = 0.0;
+  }
+
+  G4double a = 0.0, z = 0.0, amin = 0.0, amax = 0.0, amin2 = 0.0;
+  G4double amax2 = 0.0, aa = 0.0, zz = 0.0, bfis = 0.0;
+  G4double bfis0 = 0.0, ell = 0.0, el = 0.0, egs = 0.0, el80 = 0.0, el20 = 0.0;
+  G4double elmax = 0.0, sel80 = 0.0, sel20 = 0.0, x = 0.0, y = 0.0, q = 0.0, qa = 0.0, qb = 0.0;
+  G4double aj = 0.0, ak = 0.0, a1 = 0.0, a2 = 0.0;
+
+  G4int i = 0, j = 0, k = 0, m = 0;
+  G4int l = 0;
 
   G4double emncof[4][5] = {{-9.01100e+2,-1.40818e+3, 2.77000e+3,-7.06695e+2, 8.89867e+2}, 
@@ -3110 +3031 @@
 
   bfis=bfis0;
-  // assert(isnan(bfis) == false);
   
   (*sbfis)=bfis;
@@ -3155 +3075 @@
   }
 
-  // assert(isnan(elmax) == false);
   (*selmax)=elmax;
 
@@ -3164 +3083 @@
 
   x = sel20/(*selmax);
-  // assert(isnan(x) == false);
   y = sel80/(*selmax);
-  // assert(isnan(y) == false);
   
   if(el <= sel20) {
@@ -3252 +3169 @@
   // TIRAGE ALEATOIRE DANS UNE EXPONENTIELLLE : Y=EXP(-X/T)
 
-  // assert(isnan((-1*T*std::log(haz(k)))) == false);
   return (-1.0*T*std::log(haz(k)));
 }
@@ -3281 +3197 @@
 
   // ial generateur pour le cascade (et les iy pour eviter les correlations)
-  static G4int i;
+  static G4int i = 0;
   static G4int itest = 0;
   // programme principal
@@ -3288 +3204 @@
   p[pSize-1] = 8.0;
   G4double x = 0.1;
-  G4double x1;
-  G4double y;
+  G4double x1 = 0.0;
+  G4double y = 0.0;
 
   if (itest == 1) {
@@ -3312 +3228 @@
   // par regression lineaire
  fmaxhaz120:
-  standardRandom(&y, &(hazard->igraine[17]));
+    y = randomGenerator->getRandom();
+  //  standardRandom(&y, &(hazard->igraine[17]));
   i = nint(y*100);
 
@@ -3336 +3253 @@
   // Révisée pour a, z flottants 25/4/2002                           =
 
-  G4double pace2;
+  G4double pace2 = 0.0;
 
   G4int ii = idint(a+0.5);
@@ -3383 +3300 @@
   const G4int qcols = 70;
   G4double q[qrows][qcols];
+  for(G4int init_i = 0; init_i < qrows; init_i++) {
+    for(G4int init_j = 0; init_j < qcols; init_j++) {
+      q[init_i][init_j] = 0.0;
+    }
+  }
 
   G4int ix=G4int(std::floor(x+0.5));
@@ -3431 +3353 @@
   (*z_par) = z;
   (*find_par) = find;
-}
-
-
-// Fission code
-
-void G4Abla::even_odd(G4double r_origin,G4double r_even_odd,G4int &i_out)     
-{
-  // Procedure to calculate I_OUT from R_IN in a way that
-  // on the average a flat distribution in R_IN results in a
-  // fluctuating distribution in I_OUT with an even-odd effect as
-  // given by R_EVEN_ODD
-
-  //     /* ------------------------------------------------------------ */
-  //     /* EXAMPLES :                                                   */
-  //     /* ------------------------------------------------------------ */
-  //     /*    If R_EVEN_ODD = 0 :                                       */
-  //     /*           CEIL(R_IN)  ----                                   */
-  //     /*                                                              */
-  //     /*              R_IN ->                                         */
-  //     /*            (somewhere in between CEIL(R_IN) and FLOOR(R_IN)) */                                            */
-  //     /*                                                              */
-  //     /*           FLOOR(R_IN) ----       --> I_OUT                   */
-  //     /* ------------------------------------------------------------ */
-  //     /*    If R_EVEN_ODD > 0 :                                       */
-  //     /*      The interval for the above treatment is                 */
-  //     /*         larger for FLOOR(R_IN) = even and                    */
-  //     /*         smaller for FLOOR(R_IN) = odd                        */
-  //     /*    For R_EVEN_ODD < 0 : just opposite treatment              */
-  //     /* ------------------------------------------------------------ */
-
-  //     /* ------------------------------------------------------------ */
-  //     /* On input:   R_ORIGIN    nuclear charge (real number)         */
-  //     /*             R_EVEN_ODD  requested even-odd effect            */
-  //     /* Intermediate quantity: R_IN = R_ORIGIN + 0.5                 */
-  //     /* On output:  I_OUT       nuclear charge (integer)             */
-  //     /* ------------------------------------------------------------ */
-
-  //      G4double R_ORIGIN,R_IN,R_EVEN_ODD,R_REST,R_HELP;
-  G4double r_in,r_rest,r_help;
-  G4double r_floor;
-  G4double r_middle;
-  //      G4int I_OUT,N_FLOOR;
-  G4int n_floor;
-
-  r_in = r_origin + 0.5;
-  r_floor = (float)((int)(r_in));
-  if (r_even_odd < 0.001) {
-    i_out = (int)(r_floor);
-  } 
-  else {
-    r_rest = r_in - r_floor;
-    r_middle = r_floor + 0.5;
-    n_floor = (int)(r_floor);
-    if (n_floor%2 == 0) {
-      // even before modif.
-      r_help = r_middle + (r_rest - 0.5) * (1.0 - r_even_odd);
-    } 
-    else {
-      // odd before modification
-      r_help = r_middle + (r_rest - 0.5) * (1.0 + r_even_odd);
-    }
-    i_out = (int)(r_help);
-  }
-}
-
-G4double G4Abla::umass(G4double z,G4double n,G4double beta)
-{
-  // liquid-drop mass, Myers & Swiatecki, Lysekil, 1967
-  // pure liquid drop, without pairing and shell effects
-
-  // On input:    Z     nuclear charge of nucleus
-  //              N     number of neutrons in nucleus
-  //              beta  deformation of nucleus
-  // On output:   binding energy of nucleus
-
-  G4double a,umass;
-  G4double alpha;
-  G4double xcom,xvs,xe;
-  const G4double pi = 3.1416;
-     
-  a = n + z;
-  alpha = ( std::sqrt(5.0/(4.0*pi)) ) * beta;
-  // assert(isnan(alpha) == false);
-  
-  xcom = 1.0 - 1.7826 * ((a - 2.0*z)/a)*((a - 2.0*z)/a);
-  // assert(isnan(xcom) == false);
-  // factor for asymmetry dependence of surface and volume term
-  xvs = - xcom * ( 15.4941 * a - 
-                   17.9439 * std::pow(a,0.66667) * (1.0+0.4*alpha*alpha) );
-  // sum of volume and surface energy
-  xe = z*z * (0.7053/(std::pow(a,0.33333)) * (1.0-0.2*alpha*alpha) - 1.1529/a);
-  // assert(isnan(xe) == false);
-  umass = xvs + xe;
-  
-  return umass;
-}
-
-G4double G4Abla::ecoul(G4double z1,G4double n1,G4double beta1,G4double z2,G4double n2,G4double beta2,G4double d)
-{
-  // Coulomb potential between two nuclei
-  // surfaces are in a distance of d
-  // in a tip to tip configuration
-
-  // approximate formulation
-  // On input: Z1      nuclear charge of first nucleus
-  //           N1      number of neutrons in first nucleus
-  //           beta1   deformation of first nucleus
-  //           Z2      nuclear charge of second nucleus
-  //           N2      number of neutrons in second nucleus
-  //           beta2   deformation of second nucleus
-  //           d       distance of surfaces of the nuclei
-
-  //      G4double Z1,N1,beta1,Z2,N2,beta2,d,ecoul;
-  G4double ecoul;
-  G4double dtot;
-  const G4double r0 = 1.16;
-
-  dtot = r0 * ( std::pow((z1+n1),0.33333) * (1.0+(2.0/3.0)*beta1)
-                + std::pow((z2+n2),0.33333) * (1.0+(2.0/3.0)*beta2) ) + d;
-  ecoul = z1 * z2 * 1.44 / dtot;
-
-  // assert(isnan(ecoul) == false);
-  return ecoul;
-}
-
-void G4Abla::fissionDistri(G4double &a,G4double &z,G4double &e,
-                           G4double &a1,G4double &z1,G4double &e1,G4double &v1,
-                           G4double &a2,G4double &z2,G4double &e2,G4double &v2)
-{
-  //  On input: A, Z, E (mass, atomic number and exc. energy of compound nucleus
-  //                     before fission)
-  //  On output: Ai, Zi, Ei (mass, atomic number and exc. energy of fragment 1 and 2
-  //                     after fission)
-
-  //  Additionally calculated but not put in the parameter list:
-  //  Kinetic energy of prefragments EkinR1, EkinR2
-
-  //  Translation of SIMFIS18.PLI (KHS, 2.1.2001)
-
-  // This program calculates isotopic distributions of fission fragments
-  // with a semiempirical model                      
-  // Copy from SIMFIS3, KHS, 8. February 1995        
-  // Modifications made by Jose Benlliure and KHS in August 1996
-  // Energy counted from lowest barrier (J. Benlliure, KHS 1997)
-  // Some bugs corrected (J. Benlliure, KHS 1997)         
-  // Version used for thesis S. Steinhaueser (August 1997)
-  // (Curvature of LD potential increased by factor of 2!)
-
-  // Weiter veraendert mit der Absicht, eine Version zu erhalten, die
-  // derjenigen entspricht, die von J. Benlliure et al.
-  // in Nucl. Phys. A 628 (1998) 458 verwendet wurde,
-  // allerdings ohne volle Neutronenabdampfung.
-
-  // The excitation energy was calculate now for each fission channel
-  // separately. The dissipation from saddle to scission was taken from
-  // systematics, the deformation energy at scission considers the shell
-  // effects in a simplified way, and the fluctuation is included. 
-  // KHS, April 1999 
-
-  // The width in N/Z was carefully adapted to values given by Lang et al.
-
-  // The width and eventually a shift in N/Z (polarization) follows the
-  // following rules:                                        
-
-  // The line N/Z following UCD has an angle of std::atan(Zcn/Ncn)
-  // to the horizontal axis on a chart of nuclides.
-  // (For 238U the angle is 32.2 deg.)
-
-  // The following relations hold: (from Armbruster)
-  // 
-  // sigma(N) (A=const) = sigma(Z) (A=const)
-  // sigma(A) (N=const) = sigma(Z) (N=const)
-  // sigma(A) (Z=const) = sigma(N) (Z=const)
-  // 
-  // From this we get:
-  // sigma(Z) (N=const) * N = sigma(N) (Z=const) * Z
-  // sigma(A) (Z=const) = sigma(Z) (A=const) * A/Z
-  // sigma(N) (Z=const) = sigma(Z) (A=const) * A/Z
-  // Z*sigma(N) (Z=const) = N*sigma(Z) (N=const) = A*sigma(Z) (A=const)
-
-  // Excitation energy now calculated above the lowest potential point
-  // Inclusion of a distribution of excitation energies             
-
-  // Several modifications, starting from SIMFIS12: KHS November 2000 
-  // This version seems to work quite well for 238U.                  
-  // The transition from symmetric to asymmetric fission around 226Th 
-  // is reasonably well reproduced, although St. I is too strong and St. II
-  // is too weak. St. I and St. II are also weakly seen for 208Pb.
-
-  // Extensions for an event generator of fission events (21.11.2000,KHS)
-
-  // Defalt parameters (IPARS) rather carefully adjusted to
-  // pre-neutron mass distributions of Vives et al. (238U + n)
-  // Die Parameter Fgamma1 und Fgamma2 sind kleiner als die resultierenden
-  // Breiten der Massenverteilungen!!!  
-  // Fgamma1 und Fgamma2 wurden angepaￜ, so daￜ
-  // Sigma-A(ST-I) = 3.3, Sigma-A(St-II) = 5.8 (nach Vives) 
-
-  // Parameters of the model carefully adjusted by KHS (2.2.2001) to
-  // 238U + 208Pb, 1000 A MeV, Timo Enqvist et al.         
-
-
-  G4double     n;
-  G4double     nlight1,nlight2;
-  G4double     aheavy1,alight1,aheavy2,alight2;
-  G4double     eheavy1,elight1,eheavy2,elight2;
-  G4double     zheavy1_shell,zheavy2_shell;
-  G4double     zlight1,zlight2;
-  G4double     masscurv;
-  G4double     sasymm1,sasymm2,ssymm,ysum,yasymm;
-  G4double     ssymm_mode1,ssymm_mode2;
-  G4double     cz_asymm1_saddle,cz_asymm2_saddle;
-  // Curvature at saddle, modified by ld-potential
-  G4double     wzasymm1_saddle, wzasymm2_saddle, wzsymm_saddle;
-  G4double     wzasymm1_scission, wzasymm2_scission, wzsymm_scission;
-  G4double     wzasymm1,wzasymm2,wzsymm;
-  G4double     nlight1_eff, nlight2_eff;
-  G4int  imode = 0;
-  G4double     rmode;
-  G4double     z1mean = 0.0, z2mean = 0.0, z1width = 0.0, za1width = 0.0;
-  //      G4double     Z1,Z2,N1R,N2R,A1R,A2R,N1,N2,A1,A2;
-  G4double     n1r,n2r,a1r,a2r,n1,n2;
-
-  G4double     zsymm,nsymm,asymm;
-  G4double     n1mean = 0.0, n2mean, n1width;
-  G4double     dueff;
-  // effective shell effect at lowest barrier
-  G4double     eld;
-  // Excitation energy with respect to ld barrier
-  G4double     re1,re2,re3;
-  G4double     eps1,eps2;
-  G4double     n1ucd,n2ucd,z1ucd,z2ucd;
-  G4double     beta = 0.0, beta1 = 0.0, beta2 = 0.0;
-
-  G4double     dn1_pol;
-  // shift of most probable neutron number for given Z,
-  // according to polarization
-  G4int  i_help;
-
-  //   /* Parameters of the semiempirical fission model */
-  G4double a_levdens;
-  //           /* level-density parameter */
-  G4double a_levdens_light1,a_levdens_light2;
-  G4double a_levdens_heavy1,a_levdens_heavy2;
-  const G4double r_null = 1.16;
-  //          /* radius parameter */
-  G4double epsilon_1_saddle,epsilon0_1_saddle;
-  G4double epsilon_2_saddle,epsilon0_2_saddle,epsilon_symm_saddle;
-  G4double epsilon_1_scission,epsilon0_1_scission;
-  G4double epsilon_2_scission,epsilon0_2_scission;
-  G4double epsilon_symm_scission;
-  //                                   /* modified energy */
-  G4double e_eff1_saddle,e_eff2_saddle;
-  G4double epot0_mode1_saddle,epot0_mode2_saddle,epot0_symm_saddle;
-  G4double epot_mode1_saddle,epot_mode2_saddle,epot_symm_saddle;
-  G4double e_defo, e_defo1,e_defo2, e_scission = 0.0, e_asym;
-  G4double e1exc = 0.0, e2exc = 0.0;
-  G4double e1exc_sigma,e2exc_sigma;
-  G4double e1final,e2final;
-
-  const G4double r0 = 1.16;
-  G4double tker;
-  G4double ekin1,ekin2;
-  //      G4double EkinR1,EkinR2,E1,E2,V1,V2;
-  G4double ekinr1,ekinr2;
-  G4int icz,k;
-
-  //   Input parameters:
-  //OMMENT(Nuclear charge number);
-  //      G4double Z;
-  //OMMENT(Nuclear mass number);
-  //      G4double A;
-  //OMMENT(Excitation energy above fission barrier);
-  //      G4double E;
-
-  //   Model parameters:
-  //OMMENT(position of heavy peak valley 1);
-  const G4double nheavy1 = 83.0;
-  //OMMENT(position of heavy peak valley 2);
-  const G4double nheavy2 = 90.0;
-  //OMMENT(Shell effect for valley 1);
-  const G4double delta_u1_shell = -2.65;
-  //        Parameter (Delta_U1_shell = -2)
-  //OMMENT(Shell effect for valley 2);
-  const G4double delta_u2_shell = -3.8;
-  //        Parameter (Delta_U2_shell = -3.2)
-  //OMMENT(I: used shell effect);
-  G4double delta_u1;
-  //omment(I: used shell effect);
-  G4double delta_u2;
-  //OMMENT(Curvature of asymmetric valley 1);
-  const G4double cz_asymm1_shell = 0.7;
-  //OMMENT(Curvature of asymmetric valley 2);
-  const G4double cz_asymm2_shell = 0.15;
-  //OMMENT(Factor for width of distr. valley 1);
-  const G4double fwidth_asymm1 = 0.63;
-  //OMMENT(Factor for width of distr. valley 2);
-  const G4double fwidth_asymm2 = 0.97;
-  //       Parameter (CZ_asymm2_scission = 0.12)
-  //OMMENT(Parameter x: a = A/x);
-  const G4double xlevdens = 12.0;
-  //OMMENT(Factor to gamma_heavy1);
-  const G4double fgamma1 = 2.0;
-  //OMMENT(I: fading of shells (general));
-  G4double gamma;
-  //OMMENT(I: fading of shell 1);
-  G4double gamma_heavy1;
-  //OMMENT(I: fading of shell 2);
-  G4double gamma_heavy2;
-  //OMMENT(Zero-point energy at saddle);
-  const G4double e_zero_point = 0.5;
-  //OMMENT(I: friction from saddle to scission);
-  G4double e_saddle_scission;
-  //OMMENT(Friction factor);
-  const G4double friction_factor = 1.0;
-  //OMMENT(I: Internal counter for different modes); INIT(0,0,0)
-  //      Integer*4 I_MODE(3)
-  //OMMENT(I: Yield of symmetric mode);
-  G4double ysymm = 0.0;
-  //OMMENT(I: Yield of asymmetric mode 1);
-  G4double yasymm1 = 0.0;
-  //OMMENT(I: Yield of asymmetric mode 2);
-  G4double yasymm2 = 0.0;
-  //OMMENT(I: Effective position of valley 1);
-  G4double nheavy1_eff;
-  //OMMENT(I: position of heavy peak valley 1);
-  G4double zheavy1;
-  //omment(I: Effective position of valley 2);
-  G4double nheavy2_eff;
-  //OMMENT(I: position of heavy peak valley 2);
-  G4double zheavy2;
-  //omment(I: Excitation energy above saddle 1);
-  G4double eexc1_saddle;
-  //omment(I: Excitation energy above saddle 2);
-  G4double eexc2_saddle;
-  //omment(I: Excitation energy above lowest saddle);
-  G4double eexc_max;
-  //omment(I: Effective mass mode 1);
-  G4double aheavy1_mean;
-  //omment(I: Effective mass mode 2);
-  G4double aheavy2_mean;
-  //omment(I: Width of symmetric mode);
-  G4double wasymm_saddle;
-  //OMMENT(I: Width of asymmetric mode 1);
-  G4double waheavy1_saddle;
-  //OMMENT(I: Width of asymmetric mode 2);
-  G4double waheavy2_saddle;
-  //omment(I: Width of symmetric mode);
-  G4double wasymm;
-  //OMMENT(I: Width of asymmetric mode 1);
-  G4double waheavy1;
-  //OMMENT(I: Width of asymmetric mode 2);
-  G4double waheavy2;
-  //OMMENT(I: Even-odd effect in Z);
-  G4double r_e_o,r_e_o_exp;
-  //OMMENT(I: Curveture of symmetric valley);
-  G4double cz_symm;
-  //OMMENT(I: Curvature of mass distribution for fixed Z);
-  G4double cn;
-  //OMMENT(I: Curvature of Z distribution for fixed A);
-  G4double cz;
-  //OMMENT(Minimum neutron width for constant Z);
-  const G4double sigzmin = 1.16;
-  //OMMENT(Surface distance of scission configuration);
-  const G4double d = 2.0;
-
-  //   /* Charge polarisation from Wagemanns p. 397: */
-  //OMMENT(Charge polarisation standard I);
-  const G4double cpol1 = 0.65;
-  //OMMENT(Charge polarisation standard II);
-  const G4double cpol2 = 0.55;
-  //OMMENT(=1: Polarisation simult. in N and Z);
-  const G4int nzpol = 1;
-  //OMMENT(=1: test output, =0: no test output);
-  const G4int itest = 0;
-      
-  //      G4double UMASS, ECOUL, reps1, reps2, rn1_pol;
-  G4double reps1, reps2, rn1_pol;
-  //      Float_t HAZ,GAUSSHAZ;
-  G4int kkk = 0;
-  //  G4int kkk = 10; // PK
-  
-  //     I_MODE = 0;
-
-  if(itest == 1) {
-    G4cout << " cn mass " << a << G4endl;
-    G4cout << " cn charge " << z << G4endl;
-    G4cout << " cn energy " << e << G4endl;
-  }
-
-  //     /* average Z of asymmetric and symmetric components: */
-  n = a - z;  /* neutron number of the fissioning nucleus */
-
-  k = 0;
-  icz = 0;
-  if ( (std::pow(z,2)/a < 25.0) || (n < nheavy2) || (e > 500.0) ) {
-    icz = -1;
-    //          GOTO 1002;
-    goto milledeux;
-  }
-
-  nlight1 = n - nheavy1;
-  nlight2 = n - nheavy2;
-        
-  //    /* Polarisation assumed for standard I and standard II:
-  //      Z - Zucd = cpol (for A = const);
-  //      from this we get (see Armbruster)
-  //      Z - Zucd =  Acn/Ncn * cpol (for N = const)                        */
-
-  zheavy1_shell = ((nheavy1/n) * z) - ((a/n) * cpol1);
-  zheavy2_shell = ((nheavy2/n) * z) - ((a/n) * cpol2);
-
-  e_saddle_scission = 
-    (-24.0 + 0.02227 * (std::pow(z,2))/(std::pow(a,0.33333)) ) * friction_factor;
-    
-  //      /* Energy dissipated from saddle to scission                        */
-  //      /* F. Rejmund et al., Nucl. Phys. A 678 (2000) 215, fig. 4 b        */
-  //      E_saddle_scission = DMAX1(0.,E_saddle_scission);
-  if (e_saddle_scission > 0.) {
-    e_saddle_scission = e_saddle_scission;
-  }
-  else {
-    e_saddle_scission = 0.;
-  }
-  //     /* Semiempirical fission model: */
-
-  //    /* Fit to experimental result on curvature of potential at saddle */
-  //           /* reference:                                              */
-  //    /* IF Z**2/A < 33.15E0 THEN
-  //       MassCurv = 30.5438538E0 - 4.00212049E0 * Z**2/A
-  //                               + 0.11983384E0 * Z**4 / (A**2) ;
-  //     ELSE
-  //       MassCurv = 10.E0 ** (7.16993332E0 - 0.26602401E0 * Z**2/A
-  //                               + 0.00283802E0 * Z**4 / (A**2)) ;  */
-  //  /* New parametrization of T. Enqvist according to Mulgin et al. 1998 */
-  if ( (std::pow(z,2))/a < 34.0) {
-    masscurv =  std::pow( 10.0,(-1.093364 + 0.082933 * (std::pow(z,2)/a)
-                           - 0.0002602 * (std::pow(z,4)/std::pow(a,2))) );
-  } else {
-    masscurv = std::pow( 10.0,(3.053536 - 0.056477 * (std::pow(z,2)/a)
-                          + 0.0002454 * (std::pow(z,4)/std::pow(a,2))) );
-  }
-
-  cz_symm = (8.0/std::pow(z,2)) * masscurv;
-
-  if(itest == 1) {
-    G4cout << "cz_symmetry= " << cz_symm << G4endl;
-  }
-
-  if (cz_symm < 0) {
-    icz = -1;
-    //          GOTO 1002;
-    goto milledeux;
-  }
-
-  //  /* proton number in symmetric fission (centre) */
-  zsymm  = z/2.0;
-  nsymm  = n/2.0;
-  asymm = nsymm + zsymm;
-
-  zheavy1 = (cz_symm*zsymm + cz_asymm1_shell*zheavy1_shell)/(cz_symm + cz_asymm1_shell);
-  zheavy2 = (cz_symm*zsymm + cz_asymm2_shell*zheavy2_shell)/(cz_symm + cz_asymm2_shell);
-  //            /* position of valley due to influence of liquid-drop potential */
-  nheavy1_eff = (zheavy1 + (a/n * cpol1))*(n/z);
-  nheavy2_eff = (zheavy2 + (a/n * cpol2))*(n/z);
-  nlight1_eff = n - nheavy1_eff;
-  nlight2_eff = n - nheavy2_eff;
-  //  /* proton number of light fragments (centre) */
-  zlight1 = z - zheavy1;
-  //  /* proton number of light fragments (centre) */
-  zlight2 = z - zheavy2;
-  aheavy1 = nheavy1_eff + zheavy1;
-  aheavy2 = nheavy2_eff + zheavy2;
-  aheavy1_mean = aheavy1;
-  aheavy2_mean = aheavy2;
-  alight1 = nlight1_eff + zlight1;
-  alight2 = nlight2_eff + zlight2;
-
-  a_levdens = a / xlevdens;
-  a_levdens_heavy1 = aheavy1 / xlevdens;
-  a_levdens_heavy2 = aheavy2 / xlevdens;
-  a_levdens_light1 = alight1 / xlevdens;
-  a_levdens_light2 = alight2 / xlevdens;
-  gamma = a_levdens / (0.4 * (std::pow(a,1.3333)) );
-  gamma_heavy1 = ( a_levdens_heavy1 / (0.4 * (std::pow(aheavy1,1.3333)) ) ) * fgamma1;
-  gamma_heavy2 = a_levdens_heavy2 / (0.4 * (std::pow(aheavy2,1.3333)) );
-
-  cz_asymm1_saddle = cz_asymm1_shell + cz_symm;
-  cz_asymm2_saddle = cz_asymm2_shell + cz_symm;
-        
-  // Up to here: Ok! Checked CS 10/10/05           
-
-  cn = umass(zsymm,(nsymm+1.),0.0) + umass(zsymm,(nsymm-1.),0.0)
-    + 1.44 * (std::pow(zsymm,2))/
-    ( (std::pow(r_null,2)) * 
-      ( std::pow((asymm+1.0),0.33333) + std::pow((asymm-1.0),0.33333) ) *
-      ( std::pow((asymm+1.0),0.33333) + std::pow((asymm-1.0),0.33333) ) )
-    - 2.0 * umass(zsymm,nsymm,0.0)
-    - 1.44 * (std::pow(zsymm,2))/
-    ( ( 2.0 * r_null * (std::pow(asymm,0.33333)) ) * 
-      ( 2.0 * r_null * (std::pow(asymm,0.33333)) ) );
-        
-  // /* shell effect in valley of mode 1 */
-  delta_u1 = delta_u1_shell + (std::pow((zheavy1_shell-zheavy1),2))*cz_asymm1_shell;
-  // /* shell effect in valley of mode 2 */
-  delta_u2 = delta_u2_shell + (std::pow((zheavy2_shell-zheavy2),2))*cz_asymm2_shell;
-
-  //     /* liquid drop energies
-  //        at the centres of the different shell effects
-  //        with respect to liquid drop at symmetry: */
-  epot0_mode1_saddle = (std::pow((zheavy1-zsymm),2)) * cz_symm;
-  epot0_mode2_saddle = (std::pow((zheavy2-zsymm),2)) * cz_symm;
-  epot0_symm_saddle = 0.0;
-      
-  if (itest == 1) {
-    G4cout << "check zheavy1 = " << zheavy1  << G4endl;
-    G4cout << "check zheavy2 = " << zheavy2  << G4endl;
-    G4cout << "check zsymm = " << zsymm  << G4endl;
-    G4cout << "check czsymm = " << cz_symm  << G4endl;
-    G4cout << "check epot0_mode1_saddle = " << epot0_mode1_saddle  << G4endl;
-    G4cout << "check epot0_mode2_saddle = " << epot0_mode2_saddle  << G4endl;
-    G4cout << "check epot0_symm_saddle = " << epot0_symm_saddle  << G4endl;
-    G4cout << "delta_u1 = " << delta_u1 << G4endl;
-    G4cout << "delta_u2 = " << delta_u2 << G4endl;
-  }
-      
-  //     /* energies including shell effects
-  //        at the centres of the different shell effects
-  //        with respect to liquid drop at symmetry: */
-  epot_mode1_saddle = epot0_mode1_saddle + delta_u1;
-  epot_mode2_saddle = epot0_mode2_saddle + delta_u2;
-  epot_symm_saddle = epot0_symm_saddle;
-  if (itest == 1) {
-    G4cout << "check epot_mode1_saddle = " << epot_mode1_saddle  << G4endl;
-    G4cout << "check epot_mode2_saddle = " << epot_mode2_saddle  << G4endl;
-    G4cout << "check epot_symm_saddle = " << epot_symm_saddle  << G4endl;
-  }
-
-  //     /* Minimum of potential with respect to ld potential at symmetry */
-  dueff = min(epot_mode1_saddle,epot_mode2_saddle);
-  dueff = min(dueff,epot_symm_saddle);
-  dueff = dueff - epot_symm_saddle;
-
-  eld = e + dueff + e_zero_point;
-      
-  if (itest == 1) {
-    G4cout << "check dueff = " << dueff  << G4endl;
-    G4cout << "check e = " << e  << G4endl;
-    G4cout << "check e_zero_point = " << e_zero_point  << G4endl;
-    G4cout << "check eld = " << eld  << G4endl;
-  }
-  // Up to here: Ok! Checked CS 10/10/05 
-     
-  //          /* E = energy above lowest effective barrier */
-  //          /* Eld = energy above liquid-drop barrier */
-
-  //     /* Due to this treatment the energy E on input means the excitation  */
-  //     /* energy above the lowest saddle.                                   */
-
-  //  /* These energies are not used */
-  eheavy1 = e * aheavy1 / a;
-  eheavy2 = e * aheavy2 / a;
-  elight1 = e * alight1 / a;
-  elight2 = e * alight2 / a;
-
-  epsilon0_1_saddle = eld - e_zero_point - epot0_mode1_saddle;
-  //            /* excitation energy at saddle mode 1 without shell effect */
-  epsilon0_2_saddle = eld - e_zero_point - epot0_mode2_saddle;
-  //            /* excitation energy at saddle mode 2 without shell effect */
-
-  epsilon_1_saddle = eld - e_zero_point - epot_mode1_saddle;
-  //            /* excitation energy at saddle mode 1 with shell effect */
-  epsilon_2_saddle = eld - e_zero_point - epot_mode2_saddle;
-  //            /* excitation energy at saddle mode 2 with shell effect */
-  epsilon_symm_saddle = eld - e_zero_point - epot_symm_saddle;
-
-  //  /* global parameters */
-  eexc1_saddle = epsilon_1_saddle;
-  eexc2_saddle = epsilon_2_saddle;
-  eexc_max = max(eexc1_saddle,eexc2_saddle);
-  eexc_max = max(eexc_max,eld);
-       
-  //         /* EEXC_MAX is energy above the lowest saddle */
-
-
-  epsilon0_1_scission = eld + e_saddle_scission - epot0_mode1_saddle;
-  //                    /* excitation energy without shell effect */
-  epsilon0_2_scission = eld + e_saddle_scission - epot0_mode2_saddle;
-  //                    /* excitation energy without shell effect */
-
-  epsilon_1_scission = eld + e_saddle_scission - epot_mode1_saddle;
-  //                    /* excitation energy at scission */
-  epsilon_2_scission = eld+ e_saddle_scission - epot_mode2_saddle;
-  //                    /* excitation energy at scission */
-  epsilon_symm_scission = eld + e_saddle_scission - epot_symm_saddle;
-  //           /* excitation energy of symmetric fragment at scission */
-
-  //     /* Calculate widhts at the saddle:                */
-
-  e_eff1_saddle = epsilon0_1_saddle - delta_u1 * (std::exp((-epsilon_1_saddle*gamma)));
-   
-  if (e_eff1_saddle > 0.0) {
-    wzasymm1_saddle = std::sqrt( (0.5 * 
-                             (std::sqrt(1.0/a_levdens*e_eff1_saddle)) /
-                             (cz_asymm1_shell * std::exp((-epsilon_1_saddle*gamma)) + cz_symm) ) );
-  } 
-  else {
-    wzasymm1_saddle = 1.0;
-  }
-
-  e_eff2_saddle = epsilon0_2_saddle - delta_u2 * (std::exp((-epsilon_2_saddle*gamma)));
-  if (e_eff2_saddle > 0.0) {
-    wzasymm2_saddle = std::sqrt( (0.5 * 
-                             (std::sqrt(1.0/a_levdens*e_eff2_saddle)) /
-                             (cz_asymm2_shell * std::exp((-epsilon_2_saddle*gamma)) + cz_symm) ) );
-  } 
-  else {
-    wzasymm2_saddle = 1.0;
-  }
-
-  if (eld > e_zero_point) {
-    if ( (eld + epsilon_symm_saddle) < 0.0)  {
-      G4cout << "<e> eld + epsilon_symm_saddle < 0" << G4endl;
-    }
-    wzsymm_saddle = std::sqrt( (0.5 * 
-                           (std::sqrt(1.0/a_levdens*(eld+epsilon_symm_saddle))) / cz_symm ) );
-  } else {
-    wzsymm_saddle = 1.0;
-  }
-
-  if (itest == 1) {
-    G4cout << "wz1(saddle) = " << wzasymm1_saddle << G4endl;
-    G4cout << "wz2(saddle) = " << wzasymm2_saddle << G4endl;
-    G4cout << "wzsymm(saddle) = " << wzsymm_saddle << G4endl;
-  }
-     
-  //     /* Calculate widhts at the scission point: */
-  //     /* fits of ref. Beizin 1991 (Plots brought to GSI by Sergei Zhdanov) */
-
-  wzsymm_scission = wzsymm_saddle;
-
-  if (e_saddle_scission == 0.0) {
-
-    wzasymm1_scission = wzasymm1_saddle;
-    wzasymm2_scission = wzasymm2_saddle;
-
-  } 
-  else {
-
-    if (nheavy1_eff > 75.0) {
-      wzasymm1_scission = (std::sqrt(21.0)) * z/a;
-      wzasymm2_scission = (std::sqrt (max( (70.0-28.0)/3.0*(z*z/a-35.0)+28.,0.0 )) ) * z/a;
-    } 
-    else {
-      wzasymm1_scission = wzasymm1_saddle;
-      wzasymm2_scission = wzasymm2_saddle;
-    }
-
-  }
-
-  wzasymm1_scission = max(wzasymm1_scission,wzasymm1_saddle);
-  wzasymm2_scission = max(wzasymm2_scission,wzasymm2_saddle);
-
-  wzasymm1 = wzasymm1_scission * fwidth_asymm1;
-  wzasymm2 = wzasymm2_scission * fwidth_asymm2;
-  wzsymm = wzsymm_scission;
-
-  /*      if (ITEST == 1) {
-          G4cout << "WZ1(scission) = " << WZasymm1_scission << G4endl;
-          G4cout << "WZ2(scission) = " << WZasymm2_scission << G4endl;
-          G4cout << "WZsymm(scission) = " << WZsymm_scission << G4endl;
-          }
-          if (ITEST == 1) {
-          G4cout << "WZ1(scission) final= " << WZasymm1 << G4endl;
-          G4cout << "WZ2(scission) final= " << WZasymm2 << G4endl;
-          G4cout << "WZsymm(scission) final= " << WZsymm << G4endl;
-          } */
-      
-  wasymm = wzsymm * a/z;
-  waheavy1 = wzasymm1 * a/z;
-  waheavy2 = wzasymm2 * a/z;
-
-  wasymm_saddle = wzsymm_saddle * a/z;
-  waheavy1_saddle = wzasymm1_saddle * a/z;
-  waheavy2_saddle = wzasymm2_saddle * a/z;
-
-  if (itest == 1) {
-    G4cout << "wasymm = " << wzsymm << G4endl;
-    G4cout << "waheavy1 = " << waheavy1 << G4endl;
-    G4cout << "waheavy2 = " << waheavy2 << G4endl;
-  }
-  // Up to here: Ok! Checked CS 11/10/05
-            
-  if ( (epsilon0_1_saddle - delta_u1*std::exp((-epsilon_1_saddle*gamma_heavy1))) < 0.0) {
-    sasymm1 = -10.0;
-  } 
-  else {
-    sasymm1 = 2.0 * std::sqrt( a_levdens * (epsilon0_1_saddle - 
-                                       delta_u1*(std::exp((-epsilon_1_saddle*gamma_heavy1))) ) );
-  }
-
-  if ( (epsilon0_2_saddle - delta_u2*std::exp((-epsilon_2_saddle*gamma_heavy2))) < 0.0) {
-    sasymm2 = -10.0;
-  } 
-  else {
-    sasymm2 = 2.0 * std::sqrt( a_levdens * (epsilon0_2_saddle - 
-                                       delta_u2*(std::exp((-epsilon_2_saddle*gamma_heavy2))) ) );
-  }
-              
-  if (epsilon_symm_saddle > 0.0) {
-    ssymm = 2.0 * std::sqrt( a_levdens*(epsilon_symm_saddle) );
-  } 
-  else {
-    ssymm = -10.0;
-  }
-      
-  if (ssymm > -10.0) {
-    ysymm = 1.0;
-
-    if (epsilon0_1_saddle < 0.0) {
-      //  /* low energy */
-      yasymm1 = std::exp((sasymm1-ssymm)) * wzasymm1_saddle / wzsymm_saddle * 2.0;
-      //           /* factor of 2 for symmetry classes */
-    } 
-    else {
-      //        /* high energy */
-      ssymm_mode1 = 2.0 * std::sqrt( a_levdens*(epsilon0_1_saddle) );
-      yasymm1 = ( std::exp((sasymm1-ssymm)) - std::exp((ssymm_mode1 - ssymm)) )  
-        * wzasymm1_saddle / wzsymm_saddle * 2.0;
-    }
-
-    if (epsilon0_2_saddle < 0.0) {
-      //  /* low energy */
-      yasymm2 = std::exp((sasymm2-ssymm)) * wzasymm2_saddle / wzsymm_saddle * 2.0;
-      //           /* factor of 2 for symmetry classes */
-    } 
-    else {
-      //        /* high energy */
-      ssymm_mode2 = 2.0 * std::sqrt( a_levdens*(epsilon0_2_saddle) );
-      yasymm2 = ( std::exp((sasymm2-ssymm)) - std::exp((ssymm_mode2 - ssymm)) )  
-        * wzasymm2_saddle / wzsymm_saddle * 2.0;
-    }       
-    //                            /* difference in the exponent in order */
-    //                            /* to avoid numerical overflow         */
-
-  } 
-  else {
-    if ( (sasymm1 > -10.0) && (sasymm2 > -10.0) ) {
-      ysymm = 0.0;
-      yasymm1 = std::exp(sasymm1) * wzasymm1_saddle * 2.0;
-      yasymm2 = std::exp(sasymm2) * wzasymm2_saddle * 2.0;
-    }
-  }
-        
-  //  /* normalize */
-  ysum = ysymm + yasymm1 + yasymm2;
-  if (ysum > 0.0) {
-    ysymm = ysymm / ysum;
-    yasymm1 = yasymm1 / ysum;
-    yasymm2 = yasymm2 / ysum;
-    yasymm = yasymm1 + yasymm2;
-  } 
-  else {
-    ysymm = 0.0;
-    yasymm1 = 0.0;
-    yasymm2 = 0.0;
-    //        /* search minimum threshold and attribute all events to this mode */
-    if ( (epsilon_symm_saddle < epsilon_1_saddle) && (epsilon_symm_saddle < epsilon_2_saddle) ) {
-      ysymm = 1.0;
-    } 
-    else {
-      if (epsilon_1_saddle < epsilon_2_saddle) {
-        yasymm1 = 1.0;
-      } 
-      else {
-        yasymm2 = 1.0;
-      }
-    }
-  }
-
-  if (itest == 1) {
-    G4cout << "ysymm normalized= " << ysymm  << G4endl;
-    G4cout << "yasymm1 normalized= " << yasymm1  << G4endl;
-    G4cout << "yasymm2 normalized= " << yasymm2  << G4endl;
-  }
-  // Up to here: Ok! Ckecked CS 11/10/05      
-      
-  //      /* even-odd effect */
-  //      /* simple parametrization KHS, Nov. 2000. From Rejmund et al. */
-  if ((int)(z) % 2 == 0) {
-    r_e_o_exp = -0.017 * (e_saddle_scission + eld) * (e_saddle_scission + eld);
-    if ( r_e_o_exp < -307.0) {
-      r_e_o_exp = -307.0;
-      r_e_o = std::pow(10.0,r_e_o_exp);
-    }
-    else {
-      r_e_o = std::pow(10.0,r_e_o_exp);
-    }
-  } 
-  else {
-    r_e_o = 0.0;
-  }
-    
-  //      $LOOP;    /* event loop */
-  //     I_COUNT = I_COUNT + 1;
-
-  //     /* random decision: symmetric or asymmetric */
-  //     /* IMODE = 1 means asymmetric fission, mode 1,
-  //        IMODE = 2 means asymmetric fission, mode 2,
-  //        IMODE = 3 means symmetric  */
-  //      RMODE = dble(HAZ(k));
-  //      rmode = rnd.rndm();  
-  rmode = haz(k);
-  // Cast for test CS 11/10/05
-  //      RMODE = 0.54;    
-  //  rmode = 0.54;
-  if (rmode < yasymm1) {
-    imode = 1;
-  }
-  if ( (rmode > yasymm1) && (rmode < (yasymm1+yasymm2)) ) {
-    imode = 2;
-  }
-  if ( (rmode > yasymm1) && (rmode > (yasymm1+yasymm2)) ) {
-    imode = 3;
-  }
-
-  //     /* determine parameters of the Z distribution */
-  // force imode (for testing, PK)
-  // imode = 3;
-  if (imode == 1) {
-    z1mean = zheavy1;
-    z1width = wzasymm1;
-  }
-  if (imode == 2) {
-    z1mean = zheavy2;
-    z1width = wzasymm2;
-  }
-  if (imode == 3) {
-    z1mean = zsymm;
-    z1width = wzsymm;
-  }
-
-  if (itest == 1) {
-    G4cout << "nbre aleatoire tire " << rmode << G4endl;
-    G4cout << "fission mode " << imode << G4endl;
-    G4cout << "z1mean= " << z1mean << G4endl;
-    G4cout << "z1width= " << z1width << G4endl;
-  }
-                      
-  //     /* random decision: Z1 and Z2 at scission: */
-  z1 = 1.0;
-  z2 = 1.0;
-  while  ( (z1<5.0) || (z2<5.0) ) {
-    //         Z1 = dble(GAUSSHAZ(K,sngl(Z1mean),sngl(Z1width)));
-    //   z1 = rnd.gaus(z1mean,z1width);
-    z1 = gausshaz(k, z1mean, z1width);
-    z2 = z - z1;
-  }
-  if (itest == 1) {
-    G4cout << "ff charge sample " << G4endl;
-    G4cout << "z1 =  " << z1 << G4endl;
-    G4cout << "z2 = " << z2 << G4endl;
-  }
-
-  //     CALL EVEN_ODD(Z1,R_E_O,I_HELP);
-  //         /* Integer proton number with even-odd effect */
-  //     Z1 = REAL(I_HELP)
-  //      /* Z1 = INT(Z1+0.5E0); */
-  z2 = z - z1;
-
-  //     /* average N of both fragments: */
-  if (imode == 1) {
-    n1mean = (z1 + cpol1 * a/n) * n/z;
-  }
-  if (imode == 2) { 
-    n1mean = (z1 + cpol2 * a/n) * n/z;
-  }
-  /*       CASE(99)   ! only for testing;
-           N1UCD = Z1 * N/Z;
-           N2UCD = Z2 * N/Z;
-           re1 = UMASS(Z1,N1UCD,0.6) +;
-           &         UMASS(Z2,N2UCD,0.6) +;
-           &         ECOUL(Z1,N1UCD,0.6,Z2,N2UCD,0.6,d);
-           re2 = UMASS(Z1,N1UCD+1.,0.6) +;
-           &         UMASS(Z2,N2UCD-1.,0.6) +;
-           &         ECOUL(Z1,N1UCD+1.,0.6,Z2,N2UCD-1.,0.6,d);
-           re3 = UMASS(Z1,N1UCD+2.,0.6) +;
-           &         UMASS(Z2,N2UCD-2.,0.6) +;
-           &         ECOUL(Z1,N1UCD+2.,0.6,Z2,N2UCD-2.,0.6,d);
-           eps2 = (re1-2.0*re2+re3) / 2.0;
-           eps1 = re2 - re1 - eps2;
-           DN1_POL = - eps1 / (2.0 * eps2);
-           N1mean = N1UCD + DN1_POL; */
-  if (imode == 3) {
-    n1ucd = z1 * n/z;
-    n2ucd = z2 * n/z;
-    re1 = umass(z1,n1ucd,0.6) + umass(z2,n2ucd,0.6) + ecoul(z1,n1ucd,0.6,z2,n2ucd,0.6,d);
-    re2 = umass(z1,n1ucd+1.,0.6) + umass(z2,n2ucd-1.,0.6) + ecoul(z1,n1ucd+1.,0.6,z2,n2ucd-1.,0.6,d);
-    re3 = umass(z1,n1ucd+2.,0.6) + umass(z2,n2ucd-2.,0.6) + ecoul(z1,n1ucd+2.,0.6,z2,n2ucd-2.,0.6,d);
-    eps2 = (re1-2.0*re2+re3) / 2.0;
-    eps1 = re2 - re1 - eps2;
-    dn1_pol = - eps1 / (2.0 * eps2);
-    n1mean = n1ucd + dn1_pol;
-  }
-  // all fission modes features have been checked CS 11/10/05
-  n2mean = n - n1mean;
-  z2mean = z - z1mean;
-      
-  //   /* Excitation energies */
-  //     /* formulated in energies in close consistency with the fission model */
-
-  //     /* E_defo = UMASS(Z*0.5E0,N*0.5E0,0.6E0) -
-  //                 UMASS(Z*0.5E0,N*0.5E0,0);   */
-  //       /* calculates the deformation energy of the liquid drop for
-  //          deformation beta = 0.6 which is most probable at scission */
-
-  //    /* N1R and N2R provisionaly taken without fluctuations in
-  //       polarisation:                                              */
-  n1r = n1mean;
-  n2r = n2mean;
-  a1r = n1r + z1;
-  a2r = n2r + z2;
-
-  if (imode == 1) { /* N = 82 */;
-    //! /* Eexc at scission */
-    e_scission = max(epsilon_1_scission,1.0);
-    if (n1mean > (n * 0.5) ) {
-      //! /* 1. fragment is spherical */
-      beta1 = 0.0;
-      beta2 = 0.6;
-      e1exc = epsilon_1_scission * a1r / a;
-      e_defo = umass(z2,n2r,beta2) - umass(z2,n2r,0.0);
-      e2exc = epsilon_1_scission * a2r / a + e_defo;
-    }
-    else {
-      //!                       /* 2. fragment is spherical */
-      beta1 = 0.6;
-      beta2 = 0.0;
-      e_defo = umass(z1,n1r,beta1) - umass(z1,n1r,0.0);
-      e1exc = epsilon_1_scission * a1r / a + e_defo;
-      e2exc = epsilon_1_scission * a2r / a;
-    }
-  }
-           
-  if (imode == 2) { 
-    //! /* N appr. 86 */
-    e_scission = max(epsilon_2_scission,1.0);
-    if (n1mean >  (n * 0.5) ) { 
-      //! /* 2. fragment is spherical */
-      beta1 = (n1r - nheavy2) * 0.034 + 0.3;
-      e_defo = umass(z1,n1r,beta1) - umass(z1,n1r,0.0);
-      e1exc = epsilon_2_scission * a1r / a + e_defo;
-      beta2 = 0.6 - beta1;
-      e_defo = umass(z2,n2r,beta2) - umass(z2,n2r,0.0);
-      e2exc = epsilon_2_scission * a2r / a + e_defo;
-    }
-    else {
-      //!                      /* 1. fragment is spherical */
-      beta2 = (n2r - nheavy2) * 0.034 + 0.3;
-      e_defo = umass(z2,n2r,beta2) - umass(z2,n2r,0.0);
-      e2exc = epsilon_2_scission * a2r / a + e_defo;
-      beta1 = 0.6 - beta2;
-      e_defo = umass(z1,n1r,beta1) - umass(z1,n1r,0.0);
-      e1exc = epsilon_2_scission * a1r / a + e_defo;
-    }
-  }
-         
-  if (imode == 3) { 
-    // ! /* Symmetric fission channel */
-
-    //             /* the fit function for beta is the deformation for
-    //                optimum energy at the scission point, d = 2 */
-    //             /* beta  : deformation of symmetric fragments */
-    //             /* beta1 : deformation of first fragment */
-    //             /* beta2 : deformation of second fragment */
-    beta =  0.177963 + 0.0153241 * zsymm - 0.000162037 * zsymm*zsymm;
-    beta1 = 0.177963 + 0.0153241 * z1 - 0.000162037 * z1*z1;
-    //            beta1 = 0.6
-    e_defo1 = umass(z1,n1r,beta1) - umass(z1,n1r,0.0);
-    beta2 = 0.177963 + 0.0153241 * z2 - 0.000162037 * z2*z2;
-    //            beta2 = 0.6
-    e_defo2 = umass(z2,n2r,beta2) - umass(z2,n2r,0.0);
-    e_asym = umass(z1 , n1r, beta1) + umass(z2, n2r ,beta2)
-      + ecoul(z1,n1r,beta1,z2,n2r,beta2,2.0)
-      - 2.0 * umass(zsymm,nsymm,beta)
-      - ecoul(zsymm,nsymm,beta,zsymm,nsymm,beta,2.0);
-    //            E_asym = CZ_symm * (Z1 - Zsymm)**2
-    e_scission = max((epsilon_symm_scission - e_asym),1.0);
-    //         /*  $LIST(Z1,N1R,Z2,N2R,E_asym,E_scission); */
-    e1exc = e_scission * a1r / a + e_defo1;
-    e2exc = e_scission * a2r / a + e_defo2;
-  }
-  // Energies checked for all the modes CS 11/10/05
-          
-  //    /* random decision: N1R and N2R at scission, before evaporation: */
-  //    /* CN =       UMASS(Zsymm , Nsymm + 1.E0,0) +
-  //                UMASS(Zsymm, Nsymm - 1.E0,0)
-  //                 + 1.44E0 * (Zsymm)**2 /
-  //                 (r_null**2 * ((Asymm+1)**1/3 + (Asymm-1)**1/3)**2 )
-  //                 - 2.E0 * UMASS(Zsymm,Nsymm,0)
-  //                 - 1.44E0 * (Zsymm)**2 / (r_null * 2.E0 * (Asymm)**1/3)**2; */
-
-
-  //    /* N1width = std::sqrt(0.5E0 * std::sqrt(1.E0/A_levdens*(Eld+E_saddle_scission)) / CN); */
-  //    /* 8. 9. 1998: KHS (see also consideration in the first comment block)
-  //       sigma_N(Z=const) = A/Z  * sigma_Z(A=const)
-  //       sigma_Z(A=const) = 0.4 to 0.5  (from Lang paper Nucl Phys. A345 (1980) 34)
-  //       sigma_N(Z=const) = 0.45 * A/Z  (= 1.16 for 238U)
-  //       therefore: SIGZMIN = 1.16                                              */
-
-  if ( (imode == 1) || (imode == 2) ) {
-    cn=(umass(z1,n1mean+1.,beta1) + umass(z1,n1mean-1.,beta1)
-        + umass(z2,n2mean+1.,beta2) + umass(z2,n2mean-1.,beta2)
-        + ecoul(z1,n1mean+1.,beta1,z2,n2mean-1.,beta2,2.0)
-        + ecoul(z1,n1mean-1.,beta1,z2,n2mean+1.,beta2,2.0)
-        - 2.0 * ecoul(z1,n1mean,beta1,z2,n2mean,beta2,2.0)
-        - 2.0 * umass(z1, n1mean, beta1)
-        - 2.0 * umass(z2, n2mean, beta2) ) * 0.5;
-    //          /* Coulomb energy neglected for the moment! */
-    //           IF (E_scission.lt.0.) Then
-    //             write(6,*)'<E> E_scission < 0, MODE 1,2'
-    //           ENDIF
-    //           IF (CN.lt.0.) Then
-    //             write(6,*)'CN < 0, MODE 1,2'
-    //           ENDIF
-    n1width=std::sqrt( (0.5 * (std::sqrt(1.0/a_levdens*(e_scission)))/cn) );
-    n1width=max(n1width, sigzmin);
-
-    //          /* random decision: N1R and N2R at scission, before evaporation: */
-    n1r = 1.0;
-    n2r = 1.0;
-    while  ( (n1r<5.0) || (n2r<5.0) ) {
-      //             n1r = dble(gausshaz(k,sngl(n1mean),sngl(n1width)));
-      //                n1r = rnd.gaus(n1mean,n1width);
-      n1r = gausshaz(k, n1mean, n1width);
-      n2r = n - n1r;
-    }      
-    //           N1R = GAUSSHAZ(K,N1mean,N1width)
-    if (itest == 1) {
-      G4cout << "after neutron sample " << n1r << G4endl;
-    }      
-    n1r = (float)( (int)((n1r+0.5)) );
-    n2r = n - n1r;
- 
-    even_odd(z1,r_e_o,i_help);
-    //         /*  proton number with even-odd effect */
-    z1 = (float)(i_help);
-    z2 = z - z1;
-
-    a1r = z1 + n1r;
-    a2r = z2 + n2r;
-  }
-        
-  if (imode == 3) {
-    //!  /* When(3) */
-    if (nzpol > 0.0) {
-      //           /* We treat a simultaneous split in Z and N to determine polarisation */
-      cz = ( umass(z1-1., n1mean+1.,beta1)
-             +  umass(z2+1., n2mean-1.,beta1) 
-             +  umass(z1+1., n1mean-1.,beta2)
-             +  umass(z2 - 1., n2mean + 1.,beta2)
-             +  ecoul(z1-1.,n1mean+1.,beta1,z2+1.,n2mean-1.,beta2,2.0)
-             +  ecoul(z1+1.,n1mean-1.,beta1,z2-1.,n2mean+1.,beta2,2.0)
-             -  2.0 * ecoul(z1,n1mean,beta1,z2,n2mean,beta2,2.0)
-             -  2.0 * umass(z1, n1mean,beta1)
-             -  2.0 * umass(z2, n2mean,beta2) ) * 0.5;
-      //           IF (E_scission.lt.0.) Then
-      //             write(6,*) '<E> E_scission < 0, MODE 1,2'
-      //           ENDIF
-      //           IF (CZ.lt.0.) Then
-      //             write(6,*) 'CZ < 0, MODE 1,2'
-      //           ENDIF
-      za1width=std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*(e_scission)) / cz) );
-      za1width=std::sqrt( (max((za1width*za1width-(1.0/12.0)),0.1)) );
-      //                        /* Check the value of 0.1 ! */
-      //                        /* Shephard correction */
-      a1r = z1 + n1mean;
-      a1r = (float)((int)((a1r+0.5)));
-      a2r = a - a1r;
-      //           /* A1R and A2R are integer numbers now */
-      //        /* $LIST(A1R,A2R,ZA1WIDTH); */
-
-      n1ucd = n/a * a1r;
-      n2ucd = n/a * a2r;
-      z1ucd = z/a * a1r;
-      z2ucd = z/a * a2r;
-
-      re1 = umass(z1ucd-1.,n1ucd+1.,beta1) + umass(z2ucd+1.,n2ucd-1.,beta2)
-        + ecoul(z1ucd-1.,n1ucd+1.,beta1,z2ucd+1.,n2ucd-1.,beta2,d);
-      re2 = umass(z1ucd,n1ucd,beta1) + umass(z2ucd,n2ucd,beta2)
-        + ecoul(z1ucd,n1ucd,beta1,z2ucd,n2ucd,beta2,d);
-      re3 = umass(z1ucd+1.,n1ucd-1.,beta1) + umass(z2ucd-1.,n2ucd+1.,beta2) +
-        + ecoul(z1ucd+1.,n1ucd-1.,beta1,z2ucd-1.,n2ucd+1.,beta2,d);
-                 
-      eps2 = (re1-2.0*re2+re3) / 2.0;
-      eps1 = (re3 - re1)/2.0;
-      dn1_pol = - eps1 / (2.0 * eps2);
-      z1 = z1ucd + dn1_pol;
-      if (itest == 1) {
-        G4cout << "before proton sample " << z1 << G4endl;
-      }  
-      //           Z1 = dble(GAUSSHAZ(k,sngl(Z1),sngl(ZA1width)));
-      //           z1 = rnd.gaus(z1,za1width);
-      z1 = gausshaz(k, z1, za1width);
-      if (itest == 1) {
-        G4cout << "after proton sample " << z1 << G4endl;
-      }    
-      even_odd(z1,r_e_o,i_help);
-      //         /* proton number with even-odd effect */
-      z1 = (float)(i_help);
-      z2 = (float)((int)( (z - z1 + 0.5)) );
-
-      n1r = a1r - z1;
-      n2r = n - n1r;
-    } 
-    else {
-      //           /* First division of protons, then adjustment of neutrons */
-      cn = ( umass(z1, n1mean+1.,beta1) + umass(z1, n1mean-1., beta1)
-             + umass(z2, n2mean+1.,beta2) + umass(z2, n2mean-1., beta2)
-             + ecoul(z1,n1mean+1.,beta1,z2,n2mean-1.,beta2,2.0)
-             + ecoul(z1,n1mean-1.,beta1,z2,n2mean+1.,beta2,2.0)
-             - 2.0 * ecoul(z1,n1mean,beta1,z2,n2mean,beta2,2.0)
-             - 2.0 * umass(z1, n1mean, 0.6)
-             - 2.0 * umass(z2, n2mean, 0.6) ) * 0.5;
-      //          /* Coulomb energy neglected for the moment! */
-      //           IF (E_scission.lt.0.) Then
-      //             write(6,*) '<E> E_scission < 0, MODE 1,2'
-      //           Endif
-      //           IF (CN.lt.0.) Then
-      //             write(6,*) 'CN < 0, MODE 1,2'
-      //           Endif
-      n1width=std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*(e_scission)) / cn) );
-      n1width=max(n1width, sigzmin);
-
-      //          /* random decision: N1R and N2R at scission, before evaporation: */
-      //           N1R = dble(GAUSSHAZ(k,sngl(N1mean),sngl(N1width)));
-      //           n1r = rnd.gaus(n1mean,n1width);
-      n1r = gausshaz(k, n1mean, n1width);
-      n1r = (float)( (int)((n1r+0.5)) );
-      n2r = n - n1r;
-
-      even_odd(z1,r_e_o,i_help);
-      //         /* Integer proton number with even-odd effect */
-      z1 = (float)(i_help);
-      z2 = z - z1;
-
-      a1r = z1 + n1r;
-      a2r = z2 + n2r;
-          
-    }
-  }
-
-  if (itest == 1) {
-    G4cout << "remid imode = " << imode << G4endl;
-    G4cout << "n1width =  " << n1width << G4endl;
-    G4cout << "n1r = " << n1r << G4endl;
-    G4cout << "a1r = " << a1r << G4endl;
-    G4cout << "n2r = " << n2r << G4endl;
-    G4cout << "a2r = " << a2r << G4endl;
-  }
-  // Up to here: checked CS 11/10/05    
-        
-  //      /* Extracted from Lang et al. Nucl. Phys. A 345 (1980) 34 */
-  e1exc_sigma = 5.5;
-  e2exc_sigma = 5.5;
-
- neufcentquatrevingtsept:
-  //       E1final = dble(Gausshaz(k,sngl(E1exc),sngl(E1exc_sigma)));
-  //       E2final = dble(Gausshaz(k,sngl(E2exc),sngl(E2exc_sigma)));
-  //        e1final = rnd.gaus(e1exc,e1exc_sigma);
-  //        e2final = rnd.gaus(e2exc,e2exc_sigma);
-  e1final = gausshaz(k, e1exc, e1exc_sigma);
-  e2final = gausshaz(k, e2exc, e2exc_sigma);
-  if ( (e1final < 0.0) || (e2final < 0.0) ) goto neufcentquatrevingtsept;
-  if (itest == 1) {
-    G4cout << "sampled exc 1 " << e1final << G4endl;
-    G4cout << "sampled exc 2 " << e2final << G4endl;
-  }
-
-  //      /* OUTPUT QUANTITIES OF THE EVENT GENERATOR:        */
-
-  //      /* Quantities before neutron evaporation            */
-
-  //      /* Neutron number of prefragments: N1R and N2R      */
-  //      /* Atomic number of fragments: Z1 and Z2            */
-  //      /* Kinetic energy of fragments: EkinR1, EkinR2      *7
-
-  //      /* Quantities after neutron evaporation:            */
-
-  //      /* Neutron number of fragments: N1 and N2           */
-  //      /* Mass number of fragments: A1 and A2              */
-  //      /* Atomic number of fragments: Z1 and Z2            */
-  //      /* Number of evaporated neutrons: N1R-N1 and N2R-N2 */
-  //      /* Kinetic energy of fragments: EkinR1*A1/A1R and
-  //                                      EkinR2*A2/A2R       */
-
-  n1 = n1r;
-  n2 = n2r;
-  a1 = n1 + z1;
-  a2 = n2 + z2;
-  e1 = e1final;
-  e2 = e2final;
-
-  //       /* Pre-neutron-emission total kinetic energy: */
-  tker = (z1 * z2 * 1.44) /
-    ( r0 * std::pow(a1,0.33333) * (1.0 + 2.0/3.0 * beta1) +
-      r0 * std::pow(a2,0.33333) * (1.0 + 2.0/3.0 * beta2) + 2.0 );
-  //       /* Pre-neutron-emission kinetic energy of 1. fragment: */
-  ekinr1 = tker * a2 / a;
-  //       /* Pre-neutron-emission kinetic energy of 2. fragment: */
-  ekinr2 = tker * a1 / a;
-
-  v1 = std::sqrt( (ekinr1/a1) ) * 1.3887;
-  v2 = std::sqrt( (ekinr2/a2) ) * 1.3887;
-
-  if (itest == 1) {
-    G4cout << "ekinr1 " << ekinr1 << G4endl;
-    G4cout << "ekinr2 " << ekinr2 << G4endl;
-  }
-
- milledeux:       
-  //**************************
-  //*** only symmetric fission
-  //**************************
-  // Symmetric fission: Ok! Checked CS 10/10/05
-  if ( (icz == -1) || (a1 < 0.0) || (a2 < 0.0) ) {
-    //           IF (z.eq.92) THEN
-    //              write(6,*)'symmetric fission'
-    //              write(6,*)'Z,A,E,A1,A2,icz,Atot',Z,A,E,A1,A2,icz,Atot
-    //           END IF
-
-    if (itest == 1) {
-      G4cout << "milledeux: liquid-drop option "  << G4endl;
-    }
-
-    n = a-z;
-    //  proton number in symmetric fission (centre) *
-    zsymm  = z / 2.0;
-    nsymm  = n / 2.0;
-    asymm = nsymm + zsymm;
-
-    a_levdens = a / xlevdens;
-
-    masscurv = 2.0;
-    cz_symm = 8.0 / std::pow(z,2) * masscurv;
-
-    wzsymm = std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*e) / cz_symm) ) ;
-
-    if (itest == 1) {
-      G4cout << " symmetric high energy fission " << G4endl;
-      G4cout << "wzsymm " << wzsymm << G4endl;
-    }
-
-    z1mean = zsymm;
-    z1width = wzsymm;
-
-    // random decision: Z1 and Z2 at scission: */
-    z1 = 1.0;
-    z2 = 1.0;
-    while  ( (z1 < 5.0) || (z2 < 5.0) ) {
-      //           z1 = dble(gausshaz(kkk,sngl(z1mean),sngl(z1width)));
-      //           z1 = rnd.gaus(z1mean,z1width);
-      z1 = gausshaz(kkk, z1mean, z1width);
-      z2 = z - z1;
-    }
-
-    if (itest == 1) {
-      G4cout << " z1 " << z1 << G4endl;
-      G4cout << " z2 " << z2 << G4endl;
-    }
-    if (itest == 1) {
-      G4cout << " zsymm " << zsymm << G4endl;
-      G4cout << " nsymm " << nsymm << G4endl;
-      G4cout << " asymm " << asymm << G4endl;
-    }
-    //          CN =  UMASS(Zsymm , Nsymm + 1.E0) + UMASS(Zsymm, Nsymm - 1.E0)
-    //    #            + 1.44E0 * (Zsymm)**2 /
-    //    #            (r_null**2 * ((Asymm+1)**(1./3.) +
-    //    #            (Asymm-1)**(1./3.))**2 )
-    //    #            - 2.E0 * UMASS(Zsymm,Nsymm)
-    //    #            - 1.44E0 * (Zsymm)**2 /
-    //    #            (r_null * 2.E0 * (Asymm)**(1./3.))**2
-
-    n1ucd = z1 * n/z;
-    n2ucd = z2 * n/z;
-    re1 = umass(z1,n1ucd,0.6) + umass(z2,n2ucd,0.6) +
-      ecoul(z1,n1ucd,0.6,z2,n2ucd,0.6,2.0);
-    re2 = umass(z1,n1ucd+1.,0.6) + umass(z2,n2ucd-1.,0.6) +
-      ecoul(z1,n1ucd+1.,0.6,z2,n2ucd-1.,0.6,2.0);
-    re3 = umass(z1,n1ucd+2.,0.6) + umass(z2,n2ucd-2.,0.6) +
-      ecoul(z1,n1ucd+2.,0.6,z2,n2ucd-2.,0.6,2.0);
-    reps2 = (re1-2.0*re2+re3)/2.0;
-    reps1 = re2 - re1 -reps2;
-    rn1_pol = -reps1/(2.0*reps2);
-    n1mean = n1ucd + rn1_pol;
-    n2mean = n - n1mean;
-
-    if (itest == 1) {
-      G4cout << " n1mean " << n1mean << G4endl;
-      G4cout << " n2mean " << n2mean << G4endl;
-    }
-
-    cn = (umass(z1,n1mean+1.,0.0) + umass(z1,n1mean-1.,0.0) +
-          + umass(z2,n2mean+1.,0.0) + umass(z2,n2mean-1.,0.0)
-          - 2.0 * umass(z1,n1mean,0.0) +
-          - 2.0 * umass(z2,n2mean,0.0) ) * 0.5;
-    //      This is an approximation! Coulomb energy is neglected.
-
-    n1width = std::sqrt( (0.5 * std::sqrt(1.0/a_levdens*e) / cn) );
-
-    if (itest == 1) {
-      G4cout << " cn " << cn << G4endl;
-      G4cout << " n1width " << n1width << G4endl;
-    }
-        
-    // random decision: N1R and N2R at scission, before evaporation: */
-    //       N1R = dfloat(NINT(GAUSSHAZ(KKK,sngl(N1mean),sngl(N1width))));
-    //      n1r = (float)( (int)(rnd.gaus(n1mean,n1width)) );
-    n1r = (float)( (int)(gausshaz(k, n1mean,n1width)) );
-    n2r = n - n1r;
-    // Mass of first and second fragment */
-    a1 = z1 + n1r;
-    a2 = z2 + n2r;
-
-    e1 = e*a1/(a1+a2);
-    e2 = e - e*a1/(a1+a2);
-    if (itest == 1) {
-      G4cout << " n1r " << n1r << G4endl;
-      G4cout << " n2r " << n2r << G4endl;
-    }
-
-  }
-
-  if (itest == 1) {
-    G4cout << " a1 " << a1 << G4endl;
-    G4cout << " z1 " << z1 << G4endl;
-    G4cout << " a2 " << a2 << G4endl;
-    G4cout << " z2 " << z2 << G4endl;
-    G4cout << " e1 " << e1 << G4endl;
-    G4cout << " e2 " << e << G4endl;
-  }
-
-  //       /* Pre-neutron-emission total kinetic energy: */
-  tker = (z1 * z2 * 1.44) /
-    ( r0 * std::pow(a1,0.33333) * (1.0 + 2.0/3.0 * beta1) +
-      r0 * std::pow(a2,0.33333) * (1.0 + 2.0/3.0 * beta2) + 2.0 );
-  //       /* Pre-neutron-emission kinetic energy of 1. fragment: */
-  ekin1 = tker * a2 / a;
-  //       /* Pre-neutron-emission kinetic energy of 2. fragment: */
-  ekin2 = tker * a1 / a;
-
-  v1 = std::sqrt( (ekin1/a1) ) * 1.3887;
-  v2 = std::sqrt( (ekin2/a2) ) * 1.3887;
-
-  if (itest == 1) {
-    G4cout << " kinetic energies " << G4endl;
-    G4cout << " ekin1 " << ekin1 << G4endl;
-    G4cout << " ekin2 " << ekin2 << G4endl;
-  }
 }
 
@@ -4822 +3385 @@
 
   // C Matrice de rotation dans le labo:
+  G4double avv = 0.0, zvv = 0.0, enerj = 0.0, plab = 0.0, tetlab = 0.0, philab = 0.0;
+  G4int itypcasc = 0;
   G4double sitet = std::sqrt(std::pow(csrem[1],2)+std::pow(csrem[2],2));
-  G4double cstet, siphi, csphi;
+  G4double cstet = 0.0, siphi = 0.0, csphi = 0.0;
   G4double R[4][4];
-
+  for(G4int init_i = 0; init_i < 4; init_i++) {
+    for(G4int init_j = 0; init_j < 4; init_j++) {
+      R[init_i][init_j] = 0.0;
+    }
+  }
+  if(verboseLevel > 1)
+    G4cout <<"csrem = " << csrem[1] << ", " << csrem[2] << ", " << csrem[3] << ")" << G4endl;
   if(sitet > 1.0e-6) { //then
     cstet = csrem[3];
@@ -4861 +3432 @@
   G4double plabf[4];
   G4double bidon = 0.0;
-
+  for(G4int init_i = 0; init_i < 4; init_i++) {
+    plabi[init_i] = 0.0;
+    plabf[init_i] = 0.0;
+  }
+  ndec = 1;
   for(G4int i = ndec; i <= volant->iv; i++) { //do i=ndec,iv
     intp = i + nopart;
+    if(volant->copied[i]) continue; // Avoid double copying
+#ifdef USE_LEGACY_CODE
     varntp->ntrack = varntp->ntrack + 1;
+#endif
     if(nint(volant->acv[i]) == 0 && nint(volant->zpcv[i]) == 0) {
-      if(verboseLevel > 2) {
-        G4cout <<"Error: Particles with A = 0 Z = 0 detected! " << G4endl;
+      if(verboseLevel > -1) {
+        G4cout << __FILE__ << ":" << __LINE__ << " Error: Particles with A = 0 Z = 0 detected! " << G4endl;
       }
       continue;
@@ -4876 +3454 @@
       }
     }
+#ifdef USE_LEGACY_CODE
     varntp->avv[intp] = nint(volant->acv[i]);
     varntp->zvv[intp] = nint(volant->zpcv[i]);
     varntp->itypcasc[intp] = 0;    
+#else
+    avv = nint(volant->acv[i]);
+    zvv = nint(volant->zpcv[i]);
+    itypcasc = 0;    
+#endif
     // transformation de lorentz remnan --> labo:
+#ifdef USE_LEGACY_CODE
     if (varntp->avv[intp] == -1) { //then
+#else
+    if(avv == -1) {
+#endif
       masse = 138.00;   // cugnon
       // c              if (avv(intp).eq.1)  masse=938.2796     !cugnon
@@ -4887 +3475 @@
     else {
       mglms(double(volant->acv[i]),double(volant->zpcv[i]),0, &el);
-      // assert(isnan(el) == false);
       masse = volant->zpcv[i]*938.27 + (volant->acv[i] - volant->zpcv[i])*939.56 + el;
     } //end if
         
     er = std::sqrt(std::pow(volant->pcv[i],2) + std::pow(masse,2));
-    // assert(isnan(er) == false);
     plabi[1] = volant->pcv[i]*(volant->xcv[i]);
     plabi[2] = volant->pcv[i]*(volant->ycv[i]);
@@ -4898 +3484 @@
         
     ptrav2 = std::pow(plabi[1],2) + std::pow(plabi[2],2) + std::pow(plabi[3],2);
-    // assert(isnan(ptrav2) == false);
+#ifdef USE_LEGACY_CODE
     varntp->plab[intp] = std::sqrt(ptrav2); 
+#else
+    plab = std::sqrt(ptrav2);
+#endif
+#ifdef USE_LEGACY_CODE
+    if(std::abs(varntp->plab[intp] - 122.009) < 1.0e-3) {
+      G4cout <<__FILE__ << ":" << __LINE__ << " Error: varntp->plab["<< intp <<"] = " << varntp->plab[intp] << G4endl;
+      
+      volant->dump();
+      varntp->dump();
+#
+      G4cout <<"varntp->plab[intp] = " << varntp->plab[intp] << G4endl;
+      G4cout <<"ndec (starting index for loop) = " << ndec << G4endl;
+      G4cout <<"volant->iv (stopping index for loop) = " << volant->iv << G4endl;
+      G4cout <<"i (current position) = " << i << G4endl;
+      G4cout <<"intp (index for writing to varntp) = " << intp << G4endl;
+      //      exit(0);
+    }
+#endif
+
+#ifdef USE_LEGACY_CODE
     varntp->enerj[intp] = std::sqrt(ptrav2 + std::pow(masse,2)) - masse;
-
+#else
+    enerj = std::sqrt(ptrav2 + std::pow(masse,2)) - masse;
+#endif
     // Rotation dans le labo:
     for(G4int j = 1; j <= 3; j++) { //do j=1,3
       plabf[j] = 0.0;
       for(G4int k = 1; k <= 3; k++) { //do k=1,3
-        plabf[j] = plabf[j] + R[k][j]*plabi[k]; // :::Fixme::: (indices?)
+        //plabf[j] = plabf[j] + R[k][j]*plabi[k]; // :::Fixme::: (indices?)
+        plabf[j] = plabf[j] + R[j][k]*plabi[k]; // :::Fixme::: (indices?)
       } // end do
     }  // end do
     // C impulsions dans le nouveau systeme copiees dans /volant/
+#ifdef USE_LEGACY_CODE
     volant->pcv[i] = varntp->plab[intp];
+#else
+    volant->pcv[i] = plab;
+#endif
     ptrav2 = std::sqrt(std::pow(plabf[1],2) + std::pow(plabf[2],2) + std::pow(plabf[3],2));
     if(ptrav2 >= 1.0e-6) { //then
@@ -4923 +3536 @@
     } //endif
     // impulsions dans le nouveau systeme copiees dans /VAR_NTP/        
+#ifdef USE_LEGACY_CODE
     if(varntp->plab[intp] >= 1.0e-6) { //then
+#else 
+    if(plab >= 1.0e-6) { //then
+#endif
+#ifdef USE_LEGACY_CODE
       bidon = plabf[3]/(varntp->plab[intp]);
-      // assert(isnan(bidon) == false);
+#else
+      bidon = plabf[3]/plab;
+#endif
       if(bidon > 1.0) { 
         bidon = 1.0;
@@ -4932 +3552 @@
         bidon = -1.0;
       }
+#ifdef USE_LEGACY_CODE
       varntp->tetlab[intp] = std::acos(bidon);
       sitet = std::sin(varntp->tetlab[intp]);
@@ -4937 +3558 @@
       varntp->tetlab[intp] = varntp->tetlab[intp]*57.2957795;
       varntp->philab[intp] = varntp->philab[intp]*57.2957795;
+#else
+      tetlab = std::acos(bidon);
+      sitet = std::sin(tetlab);
+      philab = std::atan2(plabf[2],plabf[1]);        
+      tetlab = tetlab*57.2957795;
+      philab = philab*57.2957795;
+#endif
     }
     else {
+#ifdef USE_LEGACY_CODE
       varntp->tetlab[intp] = 90.0;
       varntp->philab[intp] = 0.0;
+#else
+      tetlab = 90.0;
+      philab = 0.0;
+#endif
     } // endif
+    volant->copied[i] = true;
+#ifndef USE_LEGACY_CODE
+    varntp->addParticle(avv, zvv, enerj, plab, tetlab, philab);
+#else
+    G4cout <<__FILE__ << ":" << __LINE__ << " volant -> varntp: " << " intp = " << intp << "Avv = " << varntp->avv[intp] << " Zvv = " << varntp->zvv[intp] << " Plab = " << varntp->plab[intp] << G4endl;
+    G4cout <<__FILE__ << ":" << __LINE__ << " volant index i = " << i << " varnpt index intp = " << intp << G4endl;
+#endif
   } // end do
 }
@@ -4971 +3611 @@
   G4double plabi[4];
   G4double plabf[4];
+  for(G4int init_i = 0; init_i < 4; init_i++) {
+    plabi[init_i] = 0.0;
+    plabf[init_i] = 0.0;
+  }
+
   // C ----Transformation de Lorentz Noyau fissionnant --> Remnant:     
   plabi[1] = p1*sitet*std::cos(phi1);
@@ -4980 +3625 @@
     plabf[j] = 0.0;
     for(G4int k = 1; k <= 3; k++) { //do k=1,3
-      //      plabf[j] = plabf[j] + R[j][k]*plabi[k];
-      plabf[j] = plabf[j] + R[k][j]*plabi[k];
+      plabf[j] = plabf[j] + R[j][k]*plabi[k];
+      //plabf[j] = plabf[j] + R[k][j]*plabi[k];
     } //end do
   } //end do
@@ -5029 +3674 @@
     pout[i] = 0.0;
     for(G4int j = 1; j <= 3; j++) { //do j=1,3
-      //      pout[i] = pout[i] + R[i][j]*pin[j];
-      pout[i] = pout[i] + R[j][i]*pin[j];
+      pout[i] = pout[i] + R[i][j]*pin[j];
+      //pout[i] = pout[i] + R[j][i]*pin[j];
     } // enddo
   } //enddo
@@ -5043 +3688 @@
   (*seed) = (*seed); // Avoid warning during compilation.
   // Use Geant4 G4UniformRand
-  (*rndm) = G4UniformRand();
+  (*rndm) = randomGenerator->getRandom();
 }
 
@@ -5050 +3695 @@
   const G4int pSize = 110;
   static G4double p[pSize];
-  static G4long ix,i;
-  static G4double x,y,a,haz;
+  static G4long ix = 0, i = 0;
+  static G4double x = 0.0, y = 0.0, a = 0.0, haz = 0.0;
   //  k =< -1 on initialise                                        
   //  k = -1 c'est reproductible                                   
@@ -5079 +3724 @@
     // the behavior of the routines in C++ and FORTRAN versions is no
     // longer exactly the same!
-    standardRandom(&x, &ix);
+    x = randomGenerator->getRandom();
+    //    standardRandom(&x, &ix);
     for(G4int i = 0; i < pSize; i++) { //do i=1,110                                                 
-      standardRandom(&(p[i]), &ix);
-    }
+      p[i] = randomGenerator->getRandom();
+      //      standardRandom(&(p[i]), &ix);
+    }
+    a = randomGenerator->getRandom();
     standardRandom(&a, &ix);
     k = 0;
@@ -5089 +3737 @@
   i = nint(100*a)+1;
   haz = p[i];
-  standardRandom(&a, &ix);
+  a = randomGenerator->getRandom();
+  //  standardRandom(&a, &ix);
   p[i] = a;
 
   hazard->ial = ix;
-  
   return haz;
 }
-
-
-G4double G4Abla::gausshaz(int k, double xmoy, double sig)
-{
-  // Gaussian random numbers:
-
-  //   1005       C*** TIRAGE ALEATOIRE DANS UNE GAUSSIENNE DE LARGEUR SIG ET MOYENNE XMOY
-  static G4int  iset = 0;
-  static G4double v1,v2,r,fac,gset,gausshaz;
-
-  if(iset == 0) { //then                                              
-    do {
-      v1 = 2.0*haz(k) - 1.0;
-      v2 = 2.0*haz(k) - 1.0;
-      r = std::pow(v1,2) + std::pow(v2,2);
-    } while(r >= 1);
-
-    fac = std::sqrt(-2.*std::log(r)/r);
-    // assert(isnan(fac) == false);
-    gset = v1*fac;
-    gausshaz = v2*fac*sig+xmoy;
-    iset = 1;
-  }
-  else {
-    gausshaz=gset*sig+xmoy;
-    iset=0;
-  }
-
-  return gausshaz;                                                         
-}
-
 
 // Utilities
@@ -5172 +3789 @@
 G4int G4Abla::nint(G4double number)
 {
-  G4double intpart;
-  G4double fractpart;
+  G4double intpart = 0.0;
+  G4double fractpart = 0.0;
   fractpart = std::modf(number, &intpart);
   if(number == 0) {

trunk/source/processes/hadronic/models/incl/src/G4AblaEvaporation.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4AblaEvaporation.cc,v 1.3 2007/10/11 08:20:08 gcosmo Exp $
+// $Id: G4AblaEvaporation.cc,v 1.4 2008/10/24 21:07:40 dennis Exp $
 //
 #include <numeric>
@@ -121 +121 @@
   G4double nucleusA = theNucleus.GetA();
   G4double nucleusZ = theNucleus.GetZ();
-  G4double nucleusMass = G4NucleiProperties::GetAtomicMass(nucleusA, nucleusZ);
+  G4double nucleusMass = G4NucleiProperties::GetNuclearMass(nucleusA, nucleusZ);
   G4double excitationEnergy = theNucleus.GetExcitationEnergy();
   G4double angularMomentum = 0.0; // Don't know how to get this quantity... From Geant4???
@@ -147 +147 @@
   }
 
-  //  G4double mTar  = G4NucleiProperties::GetAtomicMass(A, Z); // Mass of the target nucleus
+  //  G4double mTar  = G4NucleiProperties::GetNuclearMass(A, Z); // Mass of the target nucleus
   varntp->exini = exitationE;
 

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4Incl.cc,v 1.14 2007/12/04 09:49:49 gcosmo Exp $ 
+// $Id: G4Incl.cc,v 1.20 2008/11/06 10:11:27 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -33 +33 @@
 #include "G4Incl.hh"
 #include <iostream>
-#include <assert.h>
 #include "Randomize.hh"
+#include "G4InclRandomNumbers.hh"
+#include "G4Ranecu.hh"
 
 G4Incl::G4Incl()
@@ -47 +48 @@
   derivGausFunction = 4;
   densFunction = 5;
+
+  //  randomGenerator = new G4InclGeant4Random();
+  randomGenerator = new G4Ranecu();
 }
 
@@ -52 +56 @@
 {
   verboseLevel = 0;
-  
+
   // Set functions to be used for integration routine.  
   wsaxFunction = 0;
@@ -66 +70 @@
   saxw = aSaxw;
   ws = aWs;
+
+  //randomGenerator = new G4Ranecu();
+  randomGenerator = new G4InclGeant4Random();
 }
 
@@ -87 +94 @@
   varntp = aVarntp;
 
+  //  randomGenerator = new G4InclGeant4Random();
+  randomGenerator = new G4Ranecu();
   light_gaus_nuc = new G4LightGausNuc();
   light_nuc = new G4LightNuc();
@@ -118 +127 @@
 G4Incl::~G4Incl()
 {
+  delete randomGenerator;
   delete light_gaus_nuc;
   delete light_nuc;
@@ -321 +331 @@
   const G4double melec = 0.511;
 
-  G4double   pcorem;
-  G4double   pxrem;
-  G4double   pyrem;
-  G4double   pzrem;
+  G4double   pcorem = 0.0;
+  G4double   pxrem = 0.0;
+  G4double   pyrem = 0.0;
+  G4double   pzrem = 0.0;
 
   G4double ap = 0.0, zp = 0.0, mprojo = 0.0, pbeam = 0.0;
 
+  varntp->clear();
   if(calincl->f[6] == 3.0) { // pi+ 
     mprojo = 139.56995;
@@ -392 +403 @@
   G4int ibert = 1;
 
-  G4int nopart;
-  G4int izrem;
-  G4int iarem;
-  G4double esrem;
-  G4double erecrem;
-  G4double berem;
-  G4double garem;
-  G4double bimpac;
-  G4int jrem;
-  G4double alrem;
+  G4int nopart = 0;
+  G4int izrem = 0;
+  G4int iarem = 0;
+  G4double esrem = 0.0;
+  G4double erecrem = 0.0;
+  G4double berem = 0.0;
+  G4double garem = 0.0;
+  G4double bimpac = 0.0;
+  G4int jrem = 0;
+  G4double alrem = 0.0;
 
   /**
    * Coulomb barrier treatment.
    */
-  G4double probaTrans;
-  G4double rndm;
+  G4double probaTrans = 0.0;
+  G4double rndm = 0.0;
   if((calincl->f[6] == 1.0) || (calincl->f[6] >= 6.0)) {
     probaTrans = coulombTransm(calincl->f[2],ap,zp,calincl->f[0],calincl->f[1]);
@@ -656 +667 @@
   const G4double melec = 0.511;
 
-  G4double pcorem;
-  G4double pxrem;
-  G4double pyrem;
-  G4double pzrem;
+  G4double pcorem = 0.0;
+  G4double pxrem = 0.0;
+  G4double pyrem = 0.0;
+  G4double pzrem = 0.0;
 
   G4double ap = 0.0, zp = 0.0, mprojo = 0.0, pbeam = 0.0;
+
+  varntp->clear();
 
   // pi+
@@ -737 +750 @@
   G4int ibert = 1;
 
-  G4int nopart;
-  G4int izrem;
-  G4int iarem;
-  G4double esrem;
-  G4double erecrem;
-  G4double berem;
-  G4double garem;
-  G4double bimpac;
-  G4int jrem;
-  G4double alrem;
+  G4int nopart = 0;
+  G4int izrem = 0;
+  G4int iarem = 0;
+  G4double esrem = 0.0;
+  G4double erecrem = 0.0;
+  G4double berem = 0.0;
+  G4double garem = 0.0;
+  G4double bimpac = 0.0;
+  G4int jrem = 0;
+  G4double alrem = 0.0;
 
   // Coulomb barrier
   
-  G4double probaTrans;
-  G4double rndm;
+  G4double probaTrans = 0.0;
+  G4double rndm = 0.0;
 
   if((calincl->f[6] == 1.0) || (calincl->f[6] >= 6.0)) {
@@ -765 +778 @@
   // Call the actual INCL routine:
   pnu(&ibert, &nopart,&izrem,&iarem,&esrem,&erecrem,&alrem,&berem,&garem,&bimpac,&jrem);
-
+//   nopart=1;
+//   kind[0]=1;
+//   ep[0]=799.835;
+//   alpha[0]=0.08716;
+//   beta[0]=0.;
+//   gam[0]=0.99619;
+//   izrem=82;
+//   iarem=208;
+//   esrem=200.;
+//   erecrem=0.18870;
+//   alrem=-0.47101;
+//   berem=0.;
+//   garem=0.88213;
+//   bimpac=2.;
   forceAbsor(&nopart, &iarem, &izrem, &esrem, &erecrem, &alrem, &berem, &garem, &jrem);
   G4double aprf = double(iarem);    // mass number of the prefragment
@@ -836 +862 @@
         }
       }
-
+      
       if(kind[j] == 3) { 
         varntp->avv[j] = -1;
@@ -938 +964 @@
     pcorem = std::sqrt(erecrem*(erecrem + 2.0 * 938.2796 * iarem));   // cugnon
     mcorem = 938.2796 * iarem;                               // cugnon
-                
+    varntp->pcorem = pcorem;
+    varntp->mcorem = mcorem;
     // Note: Il faut negliger l'energie d'excitation (ESREM) pour que le bilan 
     // d'impulsion soit correct a la sortie de la cascade.....et prendre la
@@ -945 +972 @@
     pyrem = pcorem * berem;
     pzrem = pcorem * garem;
-        
+    varntp->pxrem = pxrem;
+    varntp->pyrem = pyrem;
+    varntp->pzrem = pzrem;
     pxbil = pxbil + pxrem;
     pybil = pybil + pyrem;
@@ -963 +992 @@
 
     // on recopie le remnant dans le ntuple
-    //  varntp->ntrack = varntp->ntrack + 1;
+    // varntp->ntrack = varntp->ntrack + 1;
     varntp->massini = iarem;
     varntp->mzini = izrem;
     varntp->exini = esrem;
 
+    if(verboseLevel > 2) {
+      G4cout << __FILE__ << ":" << __LINE__ << "Dump varntp after cascade: " << G4endl;
+      varntp->dump();
+    }
     // Evaporation/fission:
     evaporationResult->ntrack = 0;
@@ -973 +1006 @@
                     erecrem, pxrem, pyrem, pzrem, eventnumber);
 
-    //    assert(evaporationResult->ntrack > 0);
-    for(G4int evaporatedParticle = 1; evaporatedParticle < evaporationResult->ntrack; evaporatedParticle++) {
+    if(verboseLevel > 2) {
+      G4cout << __FILE__ << ":" << __LINE__ << "Dump evaporationResult after evaporation: " << G4endl;
+      evaporationResult->dump();
+    }
+    for(G4int evaporatedParticle = 0; evaporatedParticle < evaporationResult->ntrack; evaporatedParticle++) {
       if(evaporationResult->avv[evaporatedParticle] == 0 && evaporationResult->zvv[evaporatedParticle] == 0) { //Fix: Skip "empty" particles with A = 0 and Z = 0
-        continue;
+        // continue;
       }
       varntp->kfis = evaporationResult->kfis;
@@ -1009 +1045 @@
     evaporationResult->ntrack = -1;
   }
+  if(verboseLevel > 2) {
+    G4cout << __FILE__ << ":" << __LINE__ << "Dump varntp after combining: " << G4endl;
+    varntp->dump();
+  }
 }
 
@@ -1025 +1065 @@
   // input: should contain a seed (ial, odd and of 5 digits) to start the work.     
 
-  G4double xrand;
-  G4double ialdep;
-  G4int imat;
-  G4int iamat, izmat;
+  G4double xrand = 0.0;
+  G4double ialdep = 0.0;
+  G4int imat = 0;
+  G4int iamat = 0, izmat = 0;
 
   // for the 19 secondary seeds of hazard:
@@ -1045 +1085 @@
   if(initRandomSeed) {
     ialdep=hazard->ial;
-    for(int i = 0; i < IGRAINESIZE; i++) {
-      for(int j = 0; j < nbtirhaz[i]; j++) {
+    for(G4int i = 0; i < IGRAINESIZE; i++) {
+      for(G4int j = 0; j < nbtirhaz[i]; j++) {
         standardRandom(&xrand,&(hazard->ial));
       }
@@ -1095 +1135 @@
 void G4Incl::initMaterial(G4int izmat, G4int iamat, G4int imat)
 {
-  G4double res_dws;
-  G4double fnor;
-
-  G4double rcour, geom;
-  G4int nbr;
-
-  G4double step, f_r;
+  G4double res_dws = 0.0;
+  G4double fnor = 0.0;
+
+  G4double rcour = 0.0, geom = 0.0;
+  G4int nbr = 0;
+
+  G4double step = 0.0, f_r = 0.0;
 
   // rms espace r, espace p, fermi momentum and energy for light gauss nuc.      
@@ -1133 +1173 @@
   }
   
-  const G4double fmp = 938.2796;
+  const G4double fmp = 938.2796; // From INCL data
 
   // parametres moyens de densite de la cible (fermi 2 parametres)
@@ -1142 +1182 @@
   }
   else if(iamat >= 19) {
-    ws->r0 = light_nuc->r[iamat];
-    ws->adif = light_nuc->a[iamat];
+    ws->r0 = light_nuc->r[iamat-1];
+    ws->adif = light_nuc->a[iamat-1];
     ws->rmaxws = ws->r0 + (ws->xfoisa)*(ws->adif);
   }
   else if(iamat >= 6) {
-    ws->r0 = light_nuc->r[iamat];
-    ws->adif = light_nuc->a[iamat];
-    ws->rmaxws = 5.5 + 0.3*(iamat-6.)/12.;
+    ws->r0 = light_nuc->r[iamat-1];
+    ws->adif = light_nuc->a[iamat-1];
+    ws->rmaxws = 5.5 + 0.3*(double(iamat)-6.)/12.;
   }
   else if(iamat >= 2) {
@@ -1172 +1212 @@
     if(iamat == 3 && izmat == 2) {
       ws->r0 = light_gaus_nuc->rms1t[7]; // Orig: rms1t(8)
-      light_gaus_nuc->pfln[7] = light_gaus_nuc->pf1t[7]*1.291;   //!sqrt(5/3)=1.291
+      light_gaus_nuc->pfln[7] = light_gaus_nuc->pf1t[7]*1.291;   //!std::sqrt(5/3)=1.291
       light_gaus_nuc->tfln[7] = std::sqrt(std::pow(light_gaus_nuc->pfln[7],2) + fmp*fmp) - fmp;
       light_gaus_nuc->vnuc[7] = light_gaus_nuc->tfln[7] + 3.86;
@@ -1181 +1221 @@
     if(iamat == 4) {
       ws->r0 = light_gaus_nuc->rms1t[8]; // Orig: rms1t(9)
-      light_gaus_nuc->pfln[8] = light_gaus_nuc->pf1t[8]*1.291;   // !sqrt(5/3)=1.291
+      light_gaus_nuc->pfln[8] = light_gaus_nuc->pf1t[8]*1.291;   // !std::sqrt(5/3)=1.291
       light_gaus_nuc->tfln[8] = std::sqrt(std::pow(light_gaus_nuc->pfln[8],2) + fmp*fmp) - fmp;
       light_gaus_nuc->vnuc[8] = light_gaus_nuc->tfln[8] + 9.43;
@@ -1194 +1234 @@
 
   // bmax for sigma geom and various projectiles (p,n,pion/d/t/he3/he4/)
-  G4int j;
+  G4int j = 0;
   for(G4int i = 0; i < MATGEOSIZE; i++) { // Orig: do i=1,6
     j = i;
@@ -1440 +1480 @@
   G4double q[qsize];
   G4double f[qsize];
+  for(G4int init_i = 0; init_i < qsize; init_i++) {
+    q[init_i] = 0.0;
+    f[init_i] = 0.0;
+  }
+
   G4double dq=0.01;
   q[0]=0.0;
@@ -1480 +1525 @@
 G4double G4Incl::integrate(G4double ami, G4double ama, G4double step, G4int functionChoice)
 {
-  G4double res;
+  G4double res = 0.0;
   G4double x1[5];
+  for(G4int init_i = 0; init_i < 5; init_i++) {
+    x1[init_i] = 0.0;
+  }
   G4double ri = ami;
   G4double ra = ama;
-  G4int nb;
+  G4int nb = 0;
   G4double acont = 1.0;
   G4double dr = step;
@@ -1538 +1586 @@
 void G4Incl::spl2ab() 
 {
-  G4int i, j, k;
+  G4int i = 0, j = 0, k = 0;
 
   for(i=0; i <= spl2->n-3; i++) {
@@ -1561 +1609 @@
 G4double G4Incl::splineab(G4double xv)
 {
-  G4double tz;
+  G4double tz = xv-spl2->x[0];
   G4int j;
-  
-  tz = xv-spl2->x[0];
 
   if(tz < 0) {
@@ -1619 +1665 @@
   G4int l = (*l_p);
 
-  G4double minus_b1, minus_b2, minus_b3;
-  //alog 
-  G4double aml1;
-  G4double aml2; 
-  G4double amlnew; 
-  G4double arg; 
-  G4double b1;
-  G4double b2; 
-  G4double b3;
-  G4double bb2; 
-  G4double be;
+  G4double minus_b1 = 0.0, minus_b2 = 0.0, minus_b3 = 0.0;
+  G4double aml1 = 0.0;
+  G4double aml2 = 0.0;
+  G4double amlnew = 0.0;
+  G4double arg = 0.0;
+  G4double b1 = 0.0;
+  G4double b2 = 0.0;
+  G4double b3 = 0.0;
+  G4double bb2 = 0.0;
+  G4double be = 0.0;
   G4double bmass[2000]; 
-  G4double bmax2; 
-  G4double c1; 
-  G4double c2; 
-  G4double cb0; 
-  G4double cchi; 
-  G4double ccr; 
-  G4double cg; 
-  G4double cif; 
-  G4double cmultn; 
-  G4double cobe; 
-  G4double coeffb0; 
-  G4double comom;
-  G4double cstet; 
-  G4double dis1; 
-  G4double dis2; 
-  G4double dis3; 
-  G4double dist; 
-  G4double eb0; 
-  G4double ecoreh5; 
-  G4double efer; 
-  G4double egs; 
-  G4double eh5; 
-  G4double eh6; 
-  G4double eij; 
-  G4double ekout; 
-  G4double elead; 
-  G4double energie_in; 
-  G4double ener_max; 
-  G4double eout; 
+  for(G4int init_i = 0; init_i < 2000; init_i++) {
+    bmass[init_i] = 0.0;
+  }
+  G4double bmax2 = 0.0;
+  G4double c1 = 0.0;
+  G4double c2 = 0.0;
+  G4double cb0 = 0.0;
+  G4double cchi = 0.0;
+  G4double ccr = 0.0;
+  G4double cg = 0.0;
+  G4double cif = 0.0;
+  G4double cmultn = 0.0;
+  G4double cobe = 0.0;
+  G4double coeffb0 = 0.0;
+  G4double comom = 0.0;
+  G4double cstet = 0.0;
+  G4double dis1 = 0.0;
+  G4double dis2 = 0.0;
+  G4double dis3 = 0.0;
+  G4double dist = 0.0;
+  G4double eb0 = 0.0;
+  G4double ecoreh5 = 0.0;
+  G4double efer = 0.0;
+  G4double egs = 0.0;
+  G4double eh5 = 0.0;
+  G4double eh6 = 0.0;
+  G4double eij = 0.0;
+  G4double ekout = 0.0;
+  G4double elead = 0.0;
+  G4double energie_in = 0.0;
+  G4double ener_max = 0.0;
+  G4double eout = 0.0;
   G4double eps_c[BL1SIZE]; 
-  G4double epsv; 
-  G4double erecg; 
-  G4double erem; 
-  G4double exi; 
-  G4double expob0; 
-  G4double factemp; 
-  G4double fffc; 
-  G4double fm; 
-  G4double g1; 
-  G4double g2; 
-  G4double ge; 
-  G4double geff; 
-  G4double gg; 
-  G4double gl1; 
-  G4double gl2; 
-  G4double gpsg; 
-  G4int i1; 
-  G4int i20; 
-  G4int ic33; 
-  G4int ich1; 
-  G4int ich2; 
-  G4int ich3; 
-  G4int ich4; 
-  G4int ichd; 
-  G4int ichpion; 
-  G4int idecf; 
-  G4int idep; 
-  G4int iej; 
-  G4int iejn; 
-  G4int iejp; 
-  G4int i_emax; 
-  G4int iflag; 
+  for(G4int init_i = 0; init_i < BL1SIZE; init_i++) {
+    eps_c[init_i] = 0.0;
+  }
+  G4double epsv = 0.0;
+  G4double erecg = 0.0;
+  G4double erem = 0.0;
+  G4double exi = 0.0;
+  G4double expob0 = 0.0;
+  G4double factemp = 0.0;
+  G4double fffc = 0.0;
+  G4double fm = 0.0;
+  G4double g1 = 0.0;
+  G4double g2 = 0.0;
+  G4double ge = 0.0;
+  G4double geff = 0.0;
+  G4double gg = 0.0;
+  G4double gl1 = 0.0;
+  G4double gl2 = 0.0;
+  G4double gpsg = 0.0;
+  G4int i1 = 0;
+  G4int i20 = 0;
+  G4int ic33 = 0;
+  G4int ich1 = 0;
+  G4int ich2 = 0;
+  G4int ich3 = 0;
+  G4int ich4 = 0;
+  G4int ichd = 0;
+  G4int ichpion = 0;
+  G4int idecf = 0;
+  G4int idep = 0;
+  G4int iej = 0;
+  G4int iejn = 0;
+  G4int iejp = 0;
+  G4int i_emax = 0;
+  G4int iflag = 0;
   G4int iflag20 = 0; 
   G4int iflag40 = 0; 
   G4int iflag60 = 0; 
   G4int ilm = 0; 
-  G4int imin; 
+  G4int imin = 0;
   G4int indic[3000]; 
-  G4int inrem; 
-  G4int ip; 
+  for(G4int init_i = 0; init_i < 3000; init_i++) {
+    indic[init_i] = 0;
+  }
+  G4int inrem = 0;
+  G4int ip = 0;
   G4int ipi[2000]; 
-  G4int iqe; 
-  G4int irem; 
-  G4int irst_avatar; 
+  for(G4int init_i = 0; init_i < 2000; init_i++) {
+    ipi[init_i] = 0;
+  }
+  G4int iqe = 0;
+  G4int irem = 0;
+  G4int irst_avatar = 0;
   G4int isos = 0; 
-  G4int itch; 
-  G4int iteste; 
-  G4int itt; 
-  G4int ixr1; 
-  G4int ixr2; 
-  G4int ixr3; 
+  G4int itch = 0;
+  G4int iteste = 0;
+  G4int itt = 0;
+  G4int ixr1 = 0;
+  G4int ixr2 = 0;
+  G4int ixr3 = 0;
   //  G4int k; 
-  G4int kcol; 
-  G4int kd; 
+  G4int kcol = 0;
+  G4int kd = 0;
   //  G4int klm = 0; 
 //   G4int l1; 
 //   G4int l2; 
-  G4int ldel; 
-  G4int lead; 
-  G4int led; 
-  G4int lnew; 
+  G4int ldel = 0;
+  G4int lead = 0;
+  G4int led = 0;
+  G4int lnew = 0;
   G4int lp = 0; 
-  G4int lp1; 
-  G4double mcdd; 
+  G4int lp1 = 0;
+  G4double mcdd = 0.0;
   //G4double mg; 
-  G4int mg;
-  G4double mpaul1; 
-  G4double mpaul2; 
-  G4double mrdd; 
-  G4double mrdn; 
-  G4double mrdp; 
-  G4double mrnd; 
-  G4double mrnn; 
-  G4double mrpd; 
-  G4int n20; 
-  G4int nbalttf; 
-  G4int nbquit; 
-  G4int nbtest; 
-  G4int nc[300]; 
-  G4int ncol; 
-  G4int ncol_2c; 
-  G4int next; 
-  G4int nmiss; 
-  G4int np; 
-  G4int npidir; 
-  G4int npion = 0; 
-  G4int npproj[300]; 
-  G4int npx; 
-  G4int nsum_col; 
-  G4double p1v; 
-  G4double p2v; 
-  G4double p3_c[BL1SIZE]; 
-  G4double p3v; 
-  G4double pfrem1; 
-  G4double pfrem2; 
-  G4double pfrem3; 
-  G4double pfreml; 
-  G4double pfreml2; 
-  G4double phi; 
-  G4double p_mod; 
-  G4double pot; 
-  G4double pout1; 
-  G4double pout2; 
-  G4double pout3; 
-  G4double pppp; 
-  G4double prem1; 
-  G4double prem2; 
-  G4double prem3; 
-  G4double psf; 
-  G4double pspr; 
-  G4double ptotl; 
-  G4double qdeut; 
-  G4double qqq; 
-  G4double r22; 
-  G4double rcm1; 
-  G4double rcm2; 
-  G4double rcm3; 
-  G4double rcorr; 
-  G4double rhopi; 
-  G4double rndm; 
-  G4double rr; 
-  G4double rrrr; 
-  G4double s; 
-  G4double s1t1; 
-  G4double s2t1; 
-  G4double s3t1; 
-  G4double schi; 
-  G4double sepa; 
-  G4double sif; 
-  G4double sitet; 
-  G4double sp1t1; 
-  G4double sp2t1; 
-  G4double sp3t1; 
-  G4double sq = 0.0; 
-  G4double sueps; 
-  G4double t[50]; 
-  G4double t0; 
-  G4double t1; 
-  G4double t2; 
-  G4double t3; 
-  G4double t33; 
-  G4double t4; 
-  G4double t5; 
-  G4double t6; 
-  G4double t7; 
-  G4double t8; 
-  G4double tau; 
-  G4double tbid; 
-  G4double tdel; 
-  G4double temfin; 
-  G4double tim; 
-  G4double timi; 
-  G4double tlabu; 
-  G4double tp; 
-  G4double tref; 
-  G4double tri; 
-  G4double tt31; 
-  G4double tt32; 
-  G4double tt33; 
-  G4double tt34; 
-  G4double tt35; 
-  G4double tt36; 
-  G4double tte; 
-  G4double u; 
-  G4double v; 
-  G4double var_ab; 
-  G4double x; 
-  G4double x1l1; 
-  G4double x1l2; 
+  G4int mg = 0;
+  G4double mpaul1 = 0.0;
+  G4double mpaul2 = 0.0;
+  G4double mrdd = 0.0;
+  G4double mrdn = 0.0;
+  G4double mrdp = 0.0;
+  G4double mrnd = 0.0;
+  G4double mrnn = 0.0;
+  G4double mrpd = 0.0;
+  G4int n20 = 0;
+  G4int nbalttf = 0;
+  G4int nbquit = 0;
+  G4int nbtest = 0;
+  G4int nc[300];
+  G4int npproj[300];
+  for(G4int init_i = 0; init_i < 300; init_i++) {
+    nc[init_i] = 0;
+    npproj[init_i] = 0;
+  }
+  G4int ncol = 0;
+  G4int ncol_2c = 0;
+  G4int next = 0;
+  G4int nmiss = 0;
+  G4int np = 0;
+  G4int npidir = 0;
+  G4int npion = 0;
+  G4int npx = 0;
+  G4int nsum_col = 0;
+  G4double p1v = 0.0;
+  G4double p2v = 0.0;
+  G4double p3v = 0.0;
+  G4double pfrem1 = 0.0;
+  G4double pfrem2 = 0.0;
+  G4double pfrem3 = 0.0;
+  G4double pfreml = 0.0;
+  G4double pfreml2 = 0.0;
+  G4double phi = 0.0;
+  G4double p_mod = 0.0;
+  G4double pot = 0.0;
+  G4double pout1 = 0.0;
+  G4double pout2 = 0.0;
+  G4double pout3 = 0.0;
+  G4double pppp = 0.0;
+  G4double prem1 = 0.0;
+  G4double prem2 = 0.0;
+  G4double prem3 = 0.0;
+  G4double psf = 0.0;
+  G4double pspr = 0.0;
+  G4double ptotl = 0.0;
+  G4double qdeut = 0.0;
+  G4double qqq = 0.0;
+  G4double r22 = 0.0;
+  G4double rcm1 = 0.0;
+  G4double rcm2 = 0.0;
+  G4double rcm3 = 0.0;
+  G4double rcorr = 0.0;
+  G4double rhopi = 0.0;
+  G4double rndm = 0.0;
+  G4double rr = 0.0;
+  G4double rrrr = 0.0;
+  G4double s = 0.0;
+  G4double s1t1 = 0.0;
+  G4double s2t1 = 0.0;
+  G4double s3t1 = 0.0;
+  G4double schi = 0.0;
+  G4double sepa = 0.0;
+  G4double sif = 0.0;
+  G4double sitet = 0.0;
+  G4double sp1t1 = 0.0;
+  G4double sp2t1 = 0.0;
+  G4double sp3t1 = 0.0;
+  G4double sq = 0.0;
+  G4double sueps = 0.0;
+  G4double t[50];
+  for(G4int init_i = 0; init_i < 50; init_i++) {
+    t[init_i] = 0.0;
+  }
+  G4double t0 = 0.0;
+  G4double t1 = 0.0;
+  G4double t2 = 0.0;
+  G4double t3 = 0.0;
+  G4double t33 = 0.0;
+  G4double t4 = 0.0;
+  G4double t5 = 0.0;
+  G4double t6 = 0.0;
+  G4double t7 = 0.0;
+  G4double t8 = 0.0;
+  G4double tau = 0.0;
+  G4double tbid = 0.0;
+  G4double tdel = 0.0;
+  G4double temfin = 0.0;
+  G4double tim = 0.0;
+  G4double timi = 0.0;
+  G4double tlabu = 0.0;
+  G4double tp = 0.0;
+  G4double tref = 0.0;
+  G4double tri = 0.0;
+  G4double tt31 = 0.0;
+  G4double tt32 = 0.0;
+  G4double tt33 = 0.0;
+  G4double tt34 = 0.0;
+  G4double tt35 = 0.0;
+  G4double tt36 = 0.0;
+  G4double tte = 0.0;
+  G4double u = 0.0;
+  G4double v = 0.0;
+  G4double var_ab = 0.0;
+  G4double x = 0.0;
+  G4double x1l1 = 0.0;
+  G4double x1l2 = 0.0;
   G4double x1_target = 0.0; 
   G4double x2_target = 0.0; 
   G4double x3_target = 0.0; 
-  G4double x2cour; 
-  G4double x2l1; 
-  G4double x2l2; 
-  G4double x3l1; 
-  G4double x3l2; 
-  G4double xapres; 
-  G4double xavant; 
-  G4double xbl1; 
-  G4double xbl2; 
-  G4double xc; 
-  G4double xe; 
-  G4double xga; 
-  G4double xl1; 
-  G4double xl2; 
-  G4double xl3; 
-  G4double xlab; 
-  G4double xleng; 
-  G4double xlengm; 
-  G4double xmodp; 
-  G4double xpb; 
-  G4double xq; 
-  G4double xr1; 
-  G4double xr2; 
-  G4double xr3; 
-  G4double xr4; 
-  G4double xr5; 
-  G4double xr6; 
-  G4double xr7; 
-  G4double xr8; 
-  G4double xv; 
-  G4double xxx; 
-  G4double xy1; 
-  G4double xy2; 
-  G4double xy3; 
-  G4double xye; 
-  G4double y;
-  G4double q1[BL1SIZE]; 
-  G4double q2[BL1SIZE]; 
-  G4double q3[BL1SIZE]; 
-  G4double q4[BL1SIZE]; 
+  G4double x2cour = 0.0;
+  G4double x2l1 = 0.0;
+  G4double x2l2 = 0.0;
+  G4double x3l1 = 0.0;
+  G4double x3l2 = 0.0;
+  G4double xapres = 0.0;
+  G4double xavant = 0.0;
+  G4double xbl1 = 0.0;
+  G4double xbl2 = 0.0;
+  G4double xc = 0.0;
+  G4double xe = 0.0;
+  G4double xga = 0.0;
+  G4double xl1 = 0.0;
+  G4double xl2 = 0.0;
+  G4double xl3 = 0.0;
+  G4double xlab = 0.0;
+  G4double xleng = 0.0;
+  G4double xlengm = 0.0;
+  G4double xmodp = 0.0;
+  G4double xpb = 0.0;
+  G4double xq = 0.0;
+  G4double xr1 = 0.0;
+  G4double xr2 = 0.0;
+  G4double xr3 = 0.0;
+  G4double xr4 = 0.0;
+  G4double xr5 = 0.0;
+  G4double xr6 = 0.0;
+  G4double xr7 = 0.0;
+  G4double xr8 = 0.0;
+  G4double xv = 0.0;
+  G4double xxx = 0.0;
+  G4double xy1 = 0.0;
+  G4double xy2 = 0.0;
+  G4double xy3 = 0.0;
+  G4double xye = 0.0;
+  G4double y = 0.0;
+  G4double p3_c[BL1SIZE];
+  G4double q1[BL1SIZE];
+  G4double q2[BL1SIZE];
+  G4double q3[BL1SIZE];
+  G4double q4[BL1SIZE];
+  G4double ym[BL1SIZE];
+  for(G4int init_i = 0; init_i < BL1SIZE; init_i++) {
+    q1[init_i] = 0.0;
+    q2[init_i] = 0.0;
+    q3[init_i] = 0.0;
+    q4[init_i] = 0.0;
+    ym[init_i] = 0.0;
+  }
   G4double y1[BL3SIZE];
-  G4double y2[BL3SIZE]; 
-  G4double y3[BL3SIZE]; 
-  //  G4double ym[2000];
-  G4double ym[BL1SIZE]; 
-  G4double z; 
-  G4double za_i; 
-  G4double zai2; 
-  G4double zshif; 
+  G4double y2[BL3SIZE];
+  G4double y3[BL3SIZE];
+  for(G4int init_i = 0; init_i < BL1SIZE; init_i++) {
+    y1[init_i] = 0.0;
+    y2[init_i] = 0.0;
+    y3[init_i] = 0.0;
+  }
+  G4double z = 0.0;
+  G4double za_i = 0.0;
+  G4double zai2 = 0.0;
+  G4double zshif = 0.0;
   G4double ztouch = 0.0; 
-  G4double ztu; 
+  G4double ztu = 0.0;
 
   // LIEGE INC-model as a subroutine 
@@ -2209 +2284 @@
   // end for logging
 
-  //  G4int jparticip[300];
   G4int jparticip[BL1SIZE];
+  for(G4int i = 0; i < BL1SIZE; i++) {
+    jparticip[i] = 0;
+  }
+
   G4double beproj = 0.;
-  bl3->ia2 = G4int(calincl->f[0]); // f(1)->f[0] and so on..., calincl added
-  G4int iz2 = G4int(calincl->f[1]);
+  bl3->ia2 = G4int(std::floor(calincl->f[0] + 0.1)); // f(1)->f[0] and so on..., calincl added
+  G4int iz2 = G4int(std::floor(calincl->f[1] + 0.1));
   G4double r02 = 1.12;
   kindstruct->kindf7 = int(std::floor(calincl->f[6] + 0.1));
@@ -2263 +2341 @@
 
   // parametres moyens de densite de la cible (fermi 2 parametres)
-  if (bl3->ia2 > 28) { //then
+  if (bl3->ia2 >= 28) { //then
     ws->r0 = (2.745e-4*bl3->ia2+1.063)*std::pow(G4double(bl3->ia2),0.33333333);
     ws->adif = 1.63e-4*bl3->ia2 + 0.510;
     ws->rmaxws = ws->r0 + ws->xfoisa*(ws->adif);
   }
-  else if(bl3->ia2 >= 0.19) { //then
-    ws->r0 = light_nuc->r[bl3->ia2];
-    ws->adif = light_nuc->a[bl3->ia2];
+  else if(bl3->ia2 >= 19) { //then
+    ws->r0 = light_nuc->r[bl3->ia2-1];
+    ws->adif = light_nuc->a[bl3->ia2-1];
     ws->rmaxws = ws->r0 + ws->xfoisa*(ws->adif);
   }
   else if(bl3->ia2>=6) { //then
-    ws->r0 = 1.581*(light_nuc->a[bl3->ia2]) * (2.0 + 5.0 * (light_nuc->r[bl3->ia2])/(2.0 + 3.0*(light_nuc->r[bl3->ia2])));
-    ws->adif = light_nuc->a[bl3->ia2];
-    ws->rmaxws = 5.5 + 0.3*(bl3->ia2 - 6.0)/12.0;
+    ws->r0 = 1.581*(light_nuc->a[bl3->ia2-1]) * (2.0 + 5.0 * (light_nuc->r[bl3->ia2-1]))/(2.0 + 3.0*(light_nuc->r[bl3->ia2-1]));
+    ws->adif = light_nuc->a[bl3->ia2-1];
+    ws->rmaxws = 5.5 + 0.3*(double(bl3->ia2) - 6.0)/12.0;
   }
   else if(bl3->ia2 >= 2) { // then
@@ -2321 +2399 @@
   
   G4double tnr = tlab;
-  assert((tnr*tnr + 2.0*tlab*fmpinc) >= 0);
-  assert((tnr+fmpinc) != 0);
   G4double binc = std::sqrt(tnr*tnr + 2.0*tlab*fmpinc)/(tnr+fmpinc);
-  assert((1.0 - binc*binc) > 0);
   G4double ginc=1.0/std::sqrt(1.0 - binc*binc);
   G4double pinc = fmpinc*binc*ginc;
-
-//   for(G4int i = 0; i < ia; i++) {
-//     jparticip[i]=0;
-//   }
-  for(G4int i = 0; i < 300; i++) {
-    jparticip[i] = 0;
-  }
 
   for(G4int bli = 0; bli < BL1SIZE; bli++) {
@@ -2424 +2492 @@
       else {
         // deutons
-        tlabu = tlab/bl3->ia1;
+        tlabu = tlab/double(bl3->ia1);
         if (tlabu <= 400) {
           coeffb0 = -0.0035*tlabu + 6.86;
@@ -2555 +2623 @@
       standardRandom(&u, &(hazard->igraine[10]));
       cstet = u*2 - 1;
-      assert((1.0 - std::pow(cstet,2)) >= 0);
       sitet = std::sqrt(1.0 - std::pow(cstet,2));
       standardRandom(&v, &(hazard->igraine[11]));
@@ -2575 +2642 @@
 
     bl1->eps[i] = w(bl1->p1[i],bl1->p2[i],bl1->p3[i],fmp);
-    // assert(isnan(bl1->eps[i]) == false);
-    // assert(isnan(energyTest(i)) == false);
 
     sp1t1 = sp1t1 + bl1->p1[i];
@@ -2594 +2659 @@
   bl1->p3[bl3->ia1] = -sp3t1;
   bl1->eps[bl3->ia1] = w(bl1->p1[bl3->ia1],bl1->p2[bl3->ia1],bl1->p3[bl3->ia1],fmp);
-  // assert(isnan(bl1->eps[bl3->ia1]) == false);
-  // assert(isnan(energyTest(bl3->ia1)) == false);
   
   // deutons
@@ -2631 +2694 @@
       t[5] = 6.283185*t[5];   // t(6) -> t[5]                                      
       t1 = t[1];              // t(2)->t[1]                                           
-      assert((1.0 - t1*t1) >= 0);
       t2 = std::sqrt(1.0 - t1*t1);                                       
       t3 = std::cos(t[2]);  //t(3)->t[2]                                            
       t4 = std::sin(t[2]);   //t(3)->t[2]                                                                                      
       t5 = t[4];       // t(5)->t[4]                                                 
-      assert((1.0 - t5*t5) >= 0);
       t6 = std::sqrt(1.0 - t5*t5);                                         
       t7 = std::cos(t[5]);   //t(6) -> t[5]                                             
@@ -2667 +2728 @@
       bl1->ind1[i] = 0;                                                      
       bl1->eps[i] = w(bl1->p1[i],bl1->p2[i],bl1->p3[i],fmp);                            
-      // assert(isnan(bl1->eps[i]) == false);
-      // assert(isnan(energyTest(i)) == false);
       bl9->hel[i] = 0.0;
       efer = efer + bl1->eps[i] - fmp;                                        
@@ -2680 +2739 @@
     bl3->x2[ia] = -s2t1;
     bl3->x3[ia] = -s3t1;
-    // assert(isnan(bl3->x3[ia]) == false);
     bl1->p1[ia] = -sp1t1;
     bl1->p2[ia] = -sp2t1;
     bl1->p3[ia] = -sp3t1;
 
-    assert((std::pow(bl1->p1[ia],2) + std::pow(bl1->p2[ia],2) + std::pow(bl1->p3[ia],2)) >= 0);
     p_mod = std::sqrt(std::pow(bl1->p1[ia],2) + std::pow(bl1->p2[ia],2) + std::pow(bl1->p3[ia],2));
     if(p_mod > ((bl10->pf)+0.05)) {
@@ -2692 +2749 @@
 
     bl1->eps[ia] = w(bl1->p1[ia],bl1->p2[ia],bl1->p3[ia],fmp);
-    // assert(isnan(energyTest(ia)) == false);
     
     efer = efer + bl1->eps[ia]-fmp;                                        
@@ -2720 +2776 @@
       t[5] = 6.283185*t[5]; //t(6)->t[5]
       t1 = t[1]; // t(2)->t[1]
-      assert((1.0 - t1*t1) >= 0);
       t2 = std::sqrt(1.0 - t1*t1);
       t3 = std::cos(t[2]);  //t(3)->t[2]
       t4 = std::sin(t[2]); //t(3)->t[2]
       t5 = t[4]; //t(5)->t[4]
-      assert((1.0 - t5*t5) >= 0);
       t6 = std::sqrt(1.0 - t5*t5);
       t7 = std::cos(t[5]);  //t(6)->t[5]
@@ -2753 +2807 @@
       bl1->ind1[i] = 0;
       bl1->eps[i] = w(bl1->p1[i],bl1->p2[i],bl1->p3[i],fmp);
-      // assert(isnan(energyTest(i)) == false);
 
       bl9->hel[i] = 0.0;
       efer = efer + bl1->eps[i] - fmp;
     }
-    x1_target = x1_target/bl3->ia2;
-    x2_target = x2_target/bl3->ia2;
-    x3_target = x3_target/bl3->ia2;
+    x1_target = x1_target/double(bl3->ia2);
+    x2_target = x2_target/double(bl3->ia2);
+    x3_target = x3_target/double(bl3->ia2);
   }
 
@@ -2773 +2826 @@
   }
 
-  assert(z >= 0);
   z = std::sqrt(z);
   // random azimuthal direction of the impact parameter (sept 99)
@@ -2810 +2862 @@
       for(G4int i = 1; i <= bl3->ia1; i++) {
       //      for(G4int i = 1; i < bl3->ia1; i++) {
-        // assert(isnan(ginc) == false);
-        assert(ginc != 0);
         bl3->x3[i] = bl3->x3[i]/ginc;
-        // assert(isnan(bl3->x3[i]) == false);
         zai2 = ws->rmaxws*(ws->rmaxws) - std::pow((b+bl3->x1[i]),2) - std::pow(bl3->x2[i],2);
         if (zai2 < 0.0) {
@@ -2828 +2877 @@
         ilm = i;
         xlengm = xleng;
-        // assert(isnan(ztu) == false);
         ztouch = ztu;
         //      goto pnu21;
@@ -2876 +2924 @@
   b1 = (binc - be)/(1.0 - be*binc);
   b2 = -be;
-  assert((1.0 - b1*b1) > 0);
   g1 = 1.0/std::sqrt(1.0 - b1*b1);
   g2 = 1.0;
@@ -2883 +2930 @@
   if (kindstruct->kindf7 <= 2) {
     bl1->eps[1] = g1*fmp + v0;
-    assert((std::pow(bl1->eps[1],2) - std::pow(fmp,2)) >= 0);
     bl1->p3[1] = std::sqrt(std::pow(bl1->eps[1],2) - std::pow(fmp,2));
-    // assert(isnan(energyTest(1)) == false);
   }
   else {
@@ -2900 +2945 @@
       if((energie_in) <= (bl3->ia1*fmp)) {
         for(G4int i = 1; i <= bl3->ia1; i++) {
-          bl1->eps[i] = energie_in/bl3->ia1;
+          bl1->eps[i] = energie_in/double(bl3->ia1);
           bl1->p1[i] = 0.0;
           bl1->p2[i] = 0.0;
-          bl1->p3[i] = pinc/bl3->ia1;
-          // assert(isnan(energyTest(i)) == false);
+          bl1->p3[i] = pinc/double(bl3->ia1);
         }
         goto pnu1871;
@@ -2928 +2972 @@
         bl1->eps[i] = g1*(eps_c[i] + b1*p3_c[i]);
         bl1->p3[i] = g1*(b1*tte + p3_c[i]);
-        // assert(isnan(energyTest(i)) == false);
         sueps = sueps + bl1->eps[i];
       } //enddo
 
-      assert(sueps != 0);
       cobe = (tlab + fmpinc)/sueps;
 
@@ -2947 +2989 @@
       }
       for(G4int i = 1; i <= bl3->ia1; i++) { //do i=1,bl3->ia1
-        // assert(isnan(energyTest(i)) == false);
         arg = std::pow((cobe*(bl1->eps[i])),2)-pm2;
         if (arg <= 0.) { //then ! put maximum momentum to 0. 
@@ -2998 +3039 @@
       for(G4int i = 1; i <= bl3->ia1; i++) { //do i=1,bl3->ia1
       //      for(G4int i = 1; i < bl3->ia1; i++) { //do i=1,bl3->ia1
-        // assert(isnan(energyTest(i)) == false);
         arg = std::pow((cobe*(bl1->eps[i])),2) - pm2;
         comom = std::sqrt(arg/(std::pow(bl1->eps[i],2) - pm2));
-        // assert(isnan(comom) == false);
         bl1->p1[i] = comom*(bl1->p1[i]);
         bl1->p2[i] = comom*(bl1->p2[i]);
         bl1->p3[i] = comom*(bl1->p3[i]);
         bl1->eps[i] = bl1->eps[i]*cobe;
-        // assert(isnan(energyTest(i)) == false);
         if (std::fabs(am(bl1->p1[i],bl1->p2[i],bl1->p3[i],bl1->eps[i])-fmp) > 0.01) {
           if(verboseLevel > 2) {
@@ -3029 +3067 @@
       //   bl2->k = bl2->k + 1;
       //                tref=ref(x1(i),x2(i),x3(i),p1(i),p2(i),p3(i),eps(i),r22)          p-n04740
-      // assert(isnan(energyTest(i)) == false);
       tref = ref(bl3->x1[i], bl3->x2[i], bl3->x3[i], bl1->p1[i], bl1->p2[i], bl1->p3[i], bl1->eps[i], r22); 
-      // assert(isnan(tref) == false);
       if (tref > bl4->tmax5) {
         goto pnu45;
@@ -3037 +3073 @@
 
       bl2->k = bl2->k + 1;
-      assert((bl2->k >= 0) && (bl2->k < BL2INDSIZE));
       bl2->crois[bl2->k]=tref;
       bl2->ind[bl2->k]=i;
@@ -3061 +3096 @@
         }
         eij=am(bl1->p1[i]+bl1->p1[j],bl1->p2[i]+bl1->p2[j],bl1->p3[i]+bl1->p3[j],bl1->eps[i]+bl1->eps[j]);
-        // assert(isnan(eij) == false);
         if (eij < 1925.0) {
           continue;
@@ -3067 +3101 @@
         isos=bl1->ind2[i]+bl1->ind2[j];
 
-        // assert(isnan(eij) == false);
-        // assert(isnan(isos) == false);
         if (31.0*(bl3->rab2) > totalCrossSection(eij,0,isos)) {
           continue;
@@ -3074 +3106 @@
 
         bl2->k = bl2->k + 1;
-        assert((bl2->k >= 0) && (bl2->k < BL2INDSIZE));
         if (j == 1) {
           kcol = kcol + 1;
@@ -3096 +3127 @@
       // here for incoming pions:
       for(G4int i = bl3->ia1+1; i <= ia; i++) { //do i=ia1+1,ia
-        // assert(isnan(energyTest(i)) == false);
         tref = ref(bl3->x1[i], bl3->x2[i], bl3->x3[i], bl1->p1[i], bl1->p2[i], bl1->p3[i], bl1->eps[i], r22); 
-        // assert(isnan(tref) == false);
-        // assert(isnan(tref) == false);
         if (tref < bl4->tmax5) {
           bl2->k = bl2->k + 1;
-          assert((bl2->k >= 0) && (bl2->k < BL2INDSIZE));
           bl2->crois[bl2->k] = tref;
           bl2->ind[bl2->k] = i;
@@ -3123 +3150 @@
           goto pnu36;
         }
-        // assert(isnan(energyTest(i)) == false);
         tref = ref(bl3->x1[i], bl3->x2[i], bl3->x3[i], bl1->p1[i], bl1->p2[i], bl1->p3[i], bl1->eps[i], r22);
-        // assert(isnan(tref) == false);
         if(verboseLevel > 3) {
           if(tref < 0.0) {
@@ -3139 +3164 @@
         t1 = bl3->x1[i]*(bl1->p1[i])+bl3->x2[i]*(bl1->p2[i])+bl3->x3[i]*(bl1->p3[i]);                      
         t2 = bl1->p1[i]*(bl1->p1[i])+bl1->p2[i]*(bl1->p2[i])+bl1->p3[i]*(bl1->p3[i]);               
-        assert(t2 != 0);
         t3 = t1/t2;
         t4 = bl3->x1[i]*(bl3->x1[i])+bl3->x2[i]*(bl3->x2[i])+bl3->x3[i]*(bl3->x3[i]); 
@@ -3157 +3181 @@
         }
         tref = (-1.0*t3 - std::sqrt(t5))*(bl1->eps[i]);  
-        // assert(isnan(tref) == false);
         if (tref > bl4->tmax5) {
           continue;
@@ -3164 +3187 @@
       pnu37:
         bl2->k = bl2->k + 1; 
-        assert((bl2->k >= 0) && (bl2->k < BL2INDSIZE));
         bl2->crois[bl2->k] = tref; 
         bl2->ind[bl2->k] = i; 
@@ -3174 +3196 @@
       for(G4int i = bl3->ia1+1; i <= ia; i++) { //do  39 i=ia1+1,ia
         npproj[i] = 0;
-        // assert(isnan(energyTest(i)) == false);
         tref = ref(bl3->x1[i], bl3->x2[i], bl3->x3[i], bl1->p1[i], bl1->p2[i], bl1->p3[i], bl1->eps[i], r22); // line 2609 
-        // assert(isnan(tref) == false);
         if(verboseLevel > 3) {
           if(tref < 0.0) {
@@ -3198 +3218 @@
         }
         eij=am(bl1->p1[i]+bl1->p1[ilm],bl1->p2[i]+bl1->p2[ilm],bl1->p3[i]+bl1->p3[ilm],bl1->eps[i]+bl1->eps[ilm]);
-        // assert(isnan(eij) == false);
         if (eij < 1925.0) {
           continue;
         }
         isos=bl1->ind2[i]+bl1->ind2[ilm];                                             
-        // assert(isnan(eij) == false);
-        // assert(isnan(isos) == false);
         if (31.*(bl3->rab2) > totalCrossSection(eij,0,isos)) {
           continue;
@@ -3354 +3371 @@
   }
 
-  // assert(isnan(tau) == false);
   if(verboseLevel > 3) {
     G4cout <<"next = " << next << G4endl;
@@ -3521 +3537 @@
 
   // parameters for the next colliding pair
-  // assert(isnan(energyTest(l1)) == false);
-  // assert(isnan(energyTest(l2)) == false);
   t[9] = bl1->eps[bl9->l1] + bl1->eps[bl9->l2]; //t(10)->t[9] 
   t0 = 1.0/t[9]; // t(10)->t[9]
@@ -3544 +3558 @@
   }
 
-  assert(bl1->eps[bl9->l1] != 0);
   bl1->ta = tau/bl1->eps[bl9->l1];
   x1l1 = bl3->x1[bl9->l1] + bl1->p1[bl9->l1]*(bl1->ta);
@@ -3550 +3563 @@
   x3l1 = bl3->x3[bl9->l1] + bl1->p3[bl9->l1]*(bl1->ta);
 
-  assert(bl1->eps[bl9->l2] != 0);
   bl1->ta = tau/bl1->eps[bl9->l2];
   x1l2 = bl3->x1[bl9->l2] + bl1->p1[bl9->l2]*(bl1->ta);
@@ -3572 +3584 @@
   // if (k3 == 1) go to 260
   // if (k4 == 0) go to 260
-  mg=bl1->ind1[bl9->l1]+bl1->ind1[bl9->l2];
+  if(k3 != 1 && k4 != 0) {
+    mg=bl1->ind1[bl9->l1]+bl1->ind1[bl9->l2];
+    isos=bl1->ind2[bl9->l1]+bl1->ind2[bl9->l2];
+  }
   if((k3 != 1) && (k4 != 0) && (mg == 1)) {
-    isos=bl1->ind2[bl9->l1]+bl1->ind2[bl9->l2];
     // if (mg != 1) go to 260
     ldel = bl9->l2;
@@ -3580 +3594 @@
       ldel = bl9->l1;
     }
-    // assert(isnan(energyTest(ldel)) == false);
     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));
     bl6->isa = bl1->ind2[ldel];
-    // assert(isnan(sq) == false);
-    // assert(isnan(mg) == false);
-    // assert(isnan(isos) == false);
     bmax2 = totalCrossSection(sq,mg,isos)/31.415926;
     if (k5 == 0 && mg != 0) {
@@ -3593 +3603 @@
   }
   else {
-    // assert(isnan(sq) == false);
-    // assert(isnan(mg) == false);
-    // assert(isnan(isos) == false);
     bmax2 = totalCrossSection(sq,mg,isos)/31.41592;
   }
@@ -3649 +3656 @@
   if (irst_avatar == 1) {
     for(G4int i = 1; i <= bl9->l1; i = i + bl9->l1 - 1) { // bugfix!
-      assert(bl1->eps[i] != 0);
-      // assert(isnan(energyTest(i)) == false);
       bl1->ta = tau/bl1->eps[i];                                                
       bl3->x1[i] = bl3->x1[i] + bl1->p1[i]*(bl1->ta);                                     
@@ -3665 +3670 @@
   else {
     for(G4int i = 1; i <= ia; i++) {
-      // assert(isnan(energyTest(i)) == false);
-      assert(bl1->eps[i] != 0);
       bl1->ta = tau/bl1->eps[i];
       bl3->x1[i] = bl3->x1[i] + bl1->p1[i]*(bl1->ta);
@@ -3716 +3719 @@
   ich3 = bl1->ind2[bl9->l1];
   ich4 = bl1->ind2[bl9->l2];
-  // assert(isnan(energyTest(l1)) == false);
-  // assert(isnan(energyTest(l2)) == false);
   aml1 = std::sqrt(std::pow(bl1->eps[bl9->l1],2) - std::pow(bl1->p1[bl9->l1],2) - std::pow(bl1->p2[bl9->l1],2) - std::pow(bl1->p3[bl9->l1],2));
   aml2 = std::sqrt(std::pow(bl1->eps[bl9->l2],2) - std::pow(bl1->p1[bl9->l2],2) - std::pow(bl1->p2[bl9->l2],2) - std::pow(bl1->p3[bl9->l2],2));
@@ -3767 +3768 @@
     G4cout <<"Momentum: p1 = " << bl1->p1[bl9->l2] << " p2 = " << bl1->p2[bl9->l2] << " p3 = " << bl1->p3[bl9->l2] << G4endl;
   }
-  // assert(isnan(energyTest(l1)) == false);
-  // assert(isnan(energyTest(l2)) == false);
 //   bl9->l1 = l1;
 //   bl9->l2 = l2;
@@ -3777 +3776 @@
 //   l1 = bl9->l1;
 //   l2 = bl9->l2;
-  // assert(isnan(energyTest(l1)) == false);
-  // assert(isnan(energyTest(l2)) == false);
   if(verboseLevel > 3) {
     G4cout <<"End of collis call" << G4endl;
@@ -3822 +3819 @@
   }
   // restitution de l1 et l2 si rejet de la col. par pauli:
-  // assert(isnan(energyTest(l1)) == false);
-  // assert(isnan(energyTest(l2)) == false);
   bl1->p1[bl9->l1] = t[20]; //t(21)->t[20]
   bl1->p2[bl9->l1] = t[21]; //t(22)->t[21]
@@ -3836 +3831 @@
   bl1->ind2[bl9->l1] = ich3;
   bl1->ind2[bl9->l2] = ich4;
-  // assert(isnan(energyTest(l1)) == false);
-  // assert(isnan(energyTest(l2)) == false);
 
   if (bl2->k == 0) {
@@ -3896 +3889 @@
     }
 
-    if(egs < (efer-(bl3->ia2-nbalttf)*tf)) {
+    if(egs < (efer- double(bl3->ia2-nbalttf)*tf)) {
       if(varavat->kveux == 1) {
         varavat->bloc_cdpp[iavat] = 1;
@@ -3982 +3975 @@
   bl1->p3[bl9->l1] = bl1->p3[bl9->l2];
   bl1->eps[bl9->l1] = bl1->eps[bl9->l2];
-  // assert(isnan(energyTest(l1)) == false);
 
   bl3->x1[bl9->l1] = bl3->x1[bl9->l2];
@@ -3995 +3987 @@
   bl1->p3[bl9->l2] = xr3;
   bl1->eps[bl9->l2] = xr4;
-  // assert(isnan(energyTest(l2)) == false);
   
   bl3->x1[bl9->l2] = xr5;
@@ -4043 +4034 @@
       led = 1;
     }
-    if((ich1+ich2-1) == 0) {
+    if((ich1+ich2-1) > 0) {
       mcdd = mcdd + 1;
       led = 1;
@@ -4059 +4050 @@
     aml1 = am(bl1->p1[bl9->l1],bl1->p2[bl9->l1],bl1->p3[bl9->l1],bl1->eps[bl9->l1]);
     aml2 = am(bl1->p1[bl9->l2],bl1->p2[bl9->l2],bl1->p3[bl9->l2],bl1->eps[bl9->l2]);
-    // assert(isnan(aml1) == false);
-    // assert(isnan(aml2) == false);
     
     t[36] = (aml2*(bl1->p1[bl9->l1]) - aml1*(bl1->p1[bl9->l2]))/(aml1+aml2); //t(37)->t[36]
@@ -4074 +4063 @@
     t[42] = t[42]/t[45]; //t(N)->t[N-1]
     t[43] = t[43]/t[45]; //t(N)->t[N-1]
-    t[45] = t[45]/t[46];
+    t[44] = t[44]/t[45];
     cif = (t[33]*t[36] + t[34]*t[37] + t[35]*t[38])/t[39]/t[40]; //t(N)->t[N-1]
    
@@ -4103 +4092 @@
     bl1->p3[bl9->l1] = aml1*tt36/(aml1 + aml2) + t[38]; //t(39)->t[38]
     bl1->eps[bl9->l1] = w(bl1->p1[bl9->l1],bl1->p2[bl9->l1],bl1->p3[bl9->l1],aml1);
-    // assert(isnan(energyTest(l1)) == false);
     bl1->p1[bl9->l2] = aml2*tt34/(aml1 + aml2) - t[36]; //t(37)->t[36]
     bl1->p2[bl9->l2] = aml2*tt35/(aml1 + aml2) - t[37]; //t(38)->t[37]
     bl1->p3[bl9->l2] = aml2*tt36/(aml1 + aml2) - t[38]; //t(39)->t[38]
     bl1->eps[bl9->l2] = w(bl1->p1[bl9->l2],bl1->p2[bl9->l2],bl1->p3[bl9->l2],aml2);
-    // assert(isnan(energyTest(l2)) == false);
   }
   // l-conservation
@@ -4188 +4175 @@
 
   tref=ref(bl3->x1[bl9->l1], bl3->x2[bl9->l1], bl3->x3[bl9->l1], bl1->p1[bl9->l1], bl1->p2[bl9->l1], bl1->p3[bl9->l1], bl1->eps[bl9->l1],r22); // line 3502
-  // assert(isnan(tref) == false);
   
   if(verboseLevel > 3) {
@@ -4204 +4190 @@
 
   tref=ref(bl3->x1[bl9->l2], bl3->x2[bl9->l2], bl3->x3[bl9->l2], bl1->p1[bl9->l2], bl1->p2[bl9->l2], bl1->p3[bl9->l2], bl1->eps[bl9->l2],r22); // line 3516
-  // assert(isnan(tref) == false);
   
   if(verboseLevel > 3) {
@@ -4282 +4267 @@
   // largeur variable du delta (phase space factor G4introduced 4/2001)
   amlnew = std::sqrt(std::pow(bl1->eps[lnew],2)-std::pow(bl1->p1[lnew],2)-std::pow(bl1->p2[lnew],2)-std::pow(bl1->p3[lnew],2)); 
-  // assert(isnan(amlnew) == false);
   
   geff = bl1->eps[lnew]/amlnew;
   qqq = std::sqrt((std::pow(amlnew,2) - std::pow((fmp+fmpi),2))*(std::pow(amlnew,2) - std::pow((fmp-fmpi),2)))/(2.0*amlnew);
-  // assert(isnan(qqq) == false);
   
   psf = std::pow(qqq,3)/(std::pow(qqq,3) + 5832000.0);
   tdel = -hc/(g0*psf)*std::log(rndm)*geff;                                  
-  // assert(isnan(tdel) == false);
   
   if(tdel <= bl4->tmax5) {
@@ -4317 +4299 @@
   t[33] = bl1->eps[bl9->l1]; //t(34)->t[33]
   var_ab = std::pow(bl1->eps[bl9->l1],2) - std::pow(bl1->p1[bl9->l1],2) - std::pow(bl1->p2[bl9->l1],2) - std::pow(bl1->p3[bl9->l1],2);
-  assert(var_ab > 0);
   ym[npion] = 0.0;
 
@@ -4331 +4312 @@
   //  }
   // PK
-  // assert(isnan(ym[npion]) == false);
-  assert(ym[npion] != 0);
-  // assert(isnan(pcm(ym[npion], fmp, fmpi)) == false);
   if(varavat->kveux == 1) {
     varavat->del1avat[iavat] = bl1->ind1[bl9->l1];
@@ -4344 +4322 @@
     G4cout <<"q1 = " << q1[npion] << " q2 = " << q2[npion] << " q3 = " << q3[npion] << " q4 = " << q4[npion] << G4endl;
   }
-  assert(ym[npion] != 0);
   decay2(&(bl1->p1[bl9->l1]), &(bl1->p2[bl9->l1]), &(bl1->p3[bl9->l1]), &(bl1->eps[bl9->l1]), &(q1[npion]), &(q2[npion]), &(q3[npion]),
          &(q4[npion]), &(ym[npion]), &fmp, &fmpi, &(bl9->hel[bl9->l1]));
@@ -4404 +4381 @@
     }
   }
-  if(egs >= (efer-(bl3->ia2-nbalttf)*tf)) {
+  if(egs >= (efer - double(bl3->ia2-nbalttf)*tf)) {
     goto pnu850;
   }
@@ -4732 +4709 @@
   
   var_ab = std::pow(bl1->p1[bl9->l1],2) + std::pow(bl1->p2[bl9->l1],2) + std::pow(bl1->p3[bl9->l1],2);
-  assert(var_ab > 0);
   gpsg = 0.0;
   if (var_ab > 0.0) {
@@ -4798 +4774 @@
   itch = itch - (1 + bl1->ind2[bl9->l1])/2;
   var_ab = std::pow(bl1->p1[bl9->l1],2) + std::pow(bl1->p2[bl9->l1],2) + std::pow(bl1->p3[bl9->l1],2);
-  assert(var_ab > 0);
   gpsg = 0.0;
   if(var_ab > 0.0) {
@@ -4968 +4943 @@
     npion = npion + 1;
     var_ab = std::pow(bl1->eps[i],2) - std::pow(bl1->p1[i],2) - std::pow(bl1->p2[i],2) - std::pow(bl1->p3[i],2);
-    assert(var_ab > 0);
     ym[npion] = 0.0;
 
@@ -5427 +5401 @@
   G4double e1 = (*e1_p);
 
-  G4double debugOutput;
+  G4double debugOutput = 0.0;
   debugOutput = am(p1,p2,p3,e1);
-  // assert(isnan(debugOutput) == false);
 
   G4double pout11 = (*pout11_p);
@@ -5457 +5430 @@
 
   // Variables:
-  G4double a;
-  G4double aaa;
-  G4double aac;
+  G4double a = 0.0;
+  G4double aaa = 0.0;
+  G4double aac = 0.0;
   //  G4double alog;
-  G4double alphac;
+  G4double alphac = 0.0;
   //  G4double amax1;
-  G4double apt;
-  G4double argu;
-  G4double b;
-  G4double btmax;
-  G4double cfi;
-  G4double cpt;
-  G4double ctet;
-  G4double e3;
-  G4double ecm;
+  G4double apt = 0.0;
+  G4double argu = 0.0;
+  G4double b = 0.0;
+  G4double btmax = 0.0;
+  G4double cfi = 0.0;
+  G4double cpt = 0.0;
+  G4double ctet = 0.0;
+  G4double e3 = 0.0;
+  G4double ecm = 0.0;
   G4double ex[3];
   G4double ey[3];
   G4double ez[3];
-  G4double f3;
-  G4double f3max;
-  G4double fi;
-  G4double fracpn;
-  G4double heli;
-  G4int iexpi;
-  G4int ii;
-  G4int index;
-  G4int index2;
-  G4int isi;
-  G4double pin;
-  G4double pl;
-  G4double pnorm;
+  G4double qq[3];
+  for(G4int init_i = 0; init_i < 3; init_i++) {
+    ex[init_i] = 0.0;
+    ey[init_i] = 0.0;
+    ez[init_i] = 0.0;
+    qq[init_i] = 0.0;
+  }
+  G4double f3 = 0.0;
+  G4double f3max = 0.0;
+  G4double fi = 0.0;
+  G4double fracpn = 0.0;
+  G4double heli = 0.0;
+  G4int iexpi = 0;
+  G4int ii = 0;
+  G4int index = 0;
+  G4int index2 = 0;
+  G4int isi = 0;
+  G4double pin = 0.0;
+  G4double pl = 0.0;
+  G4double pnorm = 0.0;
   //  G4double pq = 0.0;
-  G4double psq;
-  G4double qq[3];
-  G4double qq4;
-  G4double ranres;
-  G4double rndm;
-  G4double s;
-  G4double s1;
+  G4double psq = 0.0;
+  G4double qq4 = 0.0;
+  G4double ranres = 0.0;
+  G4double rndm = 0.0;
+  G4double s = 0.0;
+  G4double s1 = 0.0;
   //  G4double sel;
-  G4double sfi;
-  G4double stet;
-  G4double t;
-  G4double x;
-  G4double xkh;
-  G4double xp1;
-  G4double xp2;
-  G4double xp3;
-  G4double xx;
-  G4double y;
-  G4double yn;
-  G4double z;
-  G4double zn;
-  G4double zz;
+  G4double sfi = 0.0;
+  G4double stet = 0.0;
+  G4double t = 0.0;
+  G4double x = 0.0;
+  G4double xkh = 0.0;
+  G4double xp1 = 0.0;
+  G4double xp2 = 0.0;
+  G4double xp3 = 0.0;
+  G4double xx = 0.0;
+  G4double y = 0.0;
+  G4double yn = 0.0;
+  G4double z = 0.0;
+  G4double zn = 0.0;
+  G4double zz = 0.0;
 
   // !!!  q4 = -1*q4;
@@ -5539 +5518 @@
   np = 0;
   psq = p1*p1 + p2*p2 + p3*p3;
-  assert(psq >= 0);
   pnorm = std::sqrt(psq);
   ecm = e1 + eout1;
-  // assert(isnan(ecm) == false);
   
   if(ecm < 1925.0) {
@@ -5552 +5529 @@
   }
 
-  // assert(isnan(bl8->rathr) == false);
   if(ecm < (2065.0 + bl8->rathr)) {
     goto collis17;
@@ -5579 +5555 @@
 
  collis19:
-  // assert(isnan(bl8->ramass) == false);
   if((ecm-2170.4-bl8->ramass) <= 0) {
     goto collis17;
@@ -5606 +5581 @@
   standardRandom(&rndm, &(hazard->igraine[10]));
 
-  // assert(isnan(ecm) == false);
-  // assert(isnan(iso) == false);
   s1 = lowEnergy(ecm, 1, iso);
-  // assert(isnan(s1) == false);
   
   if(m1 != 0) {
-    assert((e1*e1 - psq) >= 0);
     bl6->xx10=std::sqrt(e1*e1-psq);
     bl6->isa=is1;
-    // assert(isnan(bl6->isa) == false);
   }
   else {
-    assert((std::pow(eout1,2)-psq) >= 0);
     bl6->xx10 = std::sqrt(std::pow(eout1,2)-psq);
     bl6->isa = is2;
-    // assert(isnan(bl6->isa) == false);
   }
 
   s = s1 + srec(ecm,bl6->xx10, iso,int(bl6->isa));
-  // assert(isnan(s) == false);
-  assert(s != 0);
   a = (s - s1)/s;
-  // assert(isnan(a) == false);
   
   if((rndm-a) <= 0) {
@@ -5648 +5613 @@
   standardRandom(&rndm, &(hazard->igraine[0]));
   s = lowEnergy(ecm,0,iso);
-  // assert(isnan(ecm) == false);
-  // assert(isnan(iso) == false);
   a = deltaProductionCrossSection(ecm,iso);
   a = s/(s+a);
@@ -5662 +5625 @@
   // j.c et al, prc56(1997)2431
  collis17:
-  assert((std::pow(ecm,2) - 4.0*xm2) >= 0);
   pl = 0.5*ecm*std::sqrt(std::pow(ecm,2) - 4.0*xm2)/xm;
   x = 0.001*pl;
-  // assert(isnan(x) == false);
   
   if (iso == 0) {
@@ -5698 +5659 @@
  collis82:
   debugOutput = am(p1,p2,p3,e1);
-  // assert(isnan(debugOutput) == false);
   
   btmax = 4.0*psq*b;
@@ -5705 +5665 @@
   ranres = rndm; 
   y = 1.0 - rndm*(1.0 - z);
-  assert(y >= 0);
   t = std::log(y)/b;
   iexpi = 0;
@@ -5733 +5692 @@
         //      y = 1.0 - ranres*(10.0 - z);
         y = 1.0 - ranres*(1.0 - z);
-        assert(y >= 0);
         t = std::log(y)/alphac;
       }
@@ -5744 +5702 @@
   }      
 
-  assert((1.0 - std::pow(ctet,2)) >= 0);
   stet = std::sqrt(1.0 - std::pow(ctet,2));
   standardRandom(&rndm, &(hazard->igraine[2]));
@@ -5754 +5711 @@
 
   debugOutput = am(p1,p2,p3,e1);
-  // assert(isnan(debugOutput) == false);
 
   if(xx >= (zz*1.0e-8)) {
-    assert(xx >= 0);
     yn=std::sqrt(xx);
     zn=yn*pnorm;
-    assert(pnorm != 0);
     ez[0] = p1/pnorm; // ez(1) -> ez[0] and so on...
     ez[1] = p2/pnorm;
     ez[2] = p3/pnorm;
-    assert(yn != 0);
     ex[0] = p2/yn;
     ex[1] = -p1/yn;
     ex[2] = 0.0;     
-    assert(zn != 0);
     ey[0] = p1*p3/zn;
     ey[1] = p2*p3/zn;
@@ -5785 +5737 @@
   pout13 = -p3;
   debugOutput = am(p1,p2,p3,e1);
-  // assert(isnan(debugOutput) == false);
 
   // backward scattering according the parametrization of ref
@@ -5804 +5755 @@
 
   debugOutput = am(p1,p2,p3,e1);
-  // assert(isnan(debugOutput) == false);
   goto exitRoutine;
 
@@ -5835 +5785 @@
   y = std::pow(ecm,2);
   q2 = (y - std::pow(1076.0,2))*(y - std::pow(800.0,2))/y/4.0;                                 
-  assert(q2 >= 0);
   q3 = std::pow((std::sqrt(q2)),3);                                                  
   f3max = q3/(q3 + std::pow(180.0,3));                                      
   y = std::pow(x,2);
   q2 = (y - std::pow(1076.0,2))*(y - std::pow(800.0,2))/y/4.0;                                 
-  assert(q2 >= 0);
   q3 = std::pow((std::sqrt(q2)),3);                                                 
   f3 = q3/(q3 + std::pow(180.0,3));
@@ -5853 +5801 @@
   pin = pnorm;
   pnorm = pcm(ecm,xm,xmdel);
-  // assert(isnan(pnorm) == false);
   if (pnorm <= 0) {
     pnorm = 0.000001;
@@ -5873 +5820 @@
 
   standardRandom(&rndm, &(hazard->igraine[5]));
-  assert((std::pow(ecm,2) - 4.0*xm2) >= 0);
   x = 0.001*0.5*ecm*std::sqrt(std::pow(ecm,2) - 4.0*xm2)/xm;
   if(x < 1.4) {
@@ -5888 +5834 @@
   }
   
-  assert((1.0 - std::pow(ctet,2)) >= 0);
   stet = std::sqrt(1.0 - std::pow(ctet,2));
   standardRandom(&rndm, &(hazard->igraine[6]));
@@ -5899 +5844 @@
   zz = p3*p3;
   if(xx >= (zz*1.0e-8)) {
-    assert(xx >= 0);
     yn = std::sqrt(xx);
     zn = yn*pin;
@@ -5922 +5866 @@
   // end of correction angular distribution of delta production
 
-  assert((xp1*xp1 + xp2*xp2 + xp3*xp3 + xm*xm) >= 0);
   e3 = std::sqrt(xp1*xp1 + xp2*xp2 + xp3*xp3 + xm*xm);
   if(k4 != 0) {
@@ -5934 +5877 @@
   qq[1] = xp2;
   qq[2] = xp3;
-  assert((xp1*xp1 + xp2*xp2 + xp3*xp3 + xmdel*xmdel) >= 0);
   qq4 = std::sqrt(xp1*xp1 + xp2*xp2 + xp3*xp3 + xmdel*xmdel);
   heli = std::pow(ctet,2);
@@ -5963 +5905 @@
   }
   debugOutput = am(p1,p2,p3,e1);
-  // assert(isnan(debugOutput) == false);
 
   if (iso == 0) {
@@ -5970 +5911 @@
   if (rndm > 0.333333) {
     debugOutput = am(p1,p2,p3,e1);
-    // assert(isnan(debugOutput) == false);
     goto exitRoutine;
   }
@@ -5980 +5920 @@
   if (index == 1) {
     debugOutput = am(p1,p2,p3,e1);
-    // assert(isnan(debugOutput) == false);
     goto exitRoutine;
   }
@@ -5990 +5929 @@
   ip = -2;
   debugOutput = am(p1,p2,p3,e1);
-  // assert(isnan(debugOutput) == false);
   goto exitRoutine;
 
@@ -5998 +5936 @@
   ip = 2;
   debugOutput = am(p1,p2,p3,e1);
-  // assert(isnan(debugOutput) == false);
   goto exitRoutine;
 
@@ -6006 +5943 @@
   pout13 = -p3;
   debugOutput = am(p1,p2,p3,e1);
-  // assert(isnan(debugOutput) == false);
   goto exitRoutine;
 
@@ -6028 +5964 @@
   }
   debugOutput = am(p1,p2,p3,e1);
-  // assert(isnan(debugOutput) == false);
 
   // symmetrization of charges in pn -> n delta
@@ -6056 +5991 @@
  collis170: 
   pnorm = pcm(ecm,xm,xm);
-  // assert(isnan(pnorm) == false);
   standardRandom(&rndm, &(hazard->igraine[11]));
   ctet = -1.0 + 2.0*rndm;
@@ -6062 +5996 @@
     ctet = sign(1.0,ctet);
   }
-  assert((1.0 - ctet*ctet) >= 0);
   stet = std::sqrt(1.0 - ctet*ctet);
   standardRandom(&rndm, &(hazard->igraine[12]));
@@ -6073 +6006 @@
   m1 = 0;
   m2 = 0;
-  assert((p1*p1 + p2*p2 + p3*p3 + xm*xm) >= 0);
   e1 = std::sqrt(p1*p1 + p2*p2 + p3*p3 + xm*xm);
   eout1 = ecm - e1;
   debugOutput = am(p1,p2,p3,e1);
-  // assert(isnan(debugOutput) == false);
 
   if (iso == 0) {
@@ -6088 +6019 @@
  exitRoutine:
   debugOutput = am(p1,p2,p3,e1);
-  // assert(isnan(debugOutput) == false);
   (*p1_p) = p1;// Was pq
   (*p2_p) = p2;
@@ -6096 +6026 @@
 
   debugOutput = am(pout11,pout12,pout13,eout1);
-  // assert(isnan(debugOutput) == false);
   (*pout11_p) = pout11;
   (*pout12_p) = pout12;
@@ -6152 +6081 @@
   G4double hel = (*hel_p);
 
-  G4double rndm;
+  G4double rndm = 0.0;
 
   G4double xe = wp;
@@ -6164 +6093 @@
   // PK
   G4double xq = pcm(xi,x1,x2);                                                  
-  // assert(isnan(xq) == false);
-  G4double ctet, stet;
-
-  G4double fi, cfi, sfi;
-  G4double sal, cal;
-  G4double t1, t2;
-  G4double w1;
-  G4double beta;
+  G4double ctet = 0.0, stet = 0.0;
+
+  G4double fi = 0.0, cfi = 0.0, sfi = 0.0;
+  G4double sal = 0.0, cal = 0.0;
+  G4double t1 = 0.0, t2 = 0.0;
+  G4double w1 = 0.0;
+  G4double beta = 0.0;
 
   if(verboseLevel > 3) {
@@ -6186 +6114 @@
       ctet = sign(1.0,ctet);
     }
-    assert((1.0 - std::pow(ctet,2)) >= 0);
     stet = std::sqrt(1.0 - std::pow(ctet,2));
     standardRandom(&rndm, &(hazard->igraine[9]));
@@ -6195 +6122 @@
   cfi = std::cos(fi);
   sfi = std::sin(fi);
-  assert((b1*b1+b2*b2+b3*b3) >= 0);
   beta = std::sqrt(b1*b1+b2*b2+b3*b3);
-  // assert(isnan(beta) == false);
-  assert(beta != 0);
-
-  assert((std::pow(b1,2) + std::pow(b2,2)) >= 0);
+
   sal = std::sqrt(std::pow(b1,2) + std::pow(b2,2))/beta;
   cal = b3/beta;
@@ -6206 +6129 @@
   if((beta >= 1.0e-10) || (sal >= 1.0e-6)) {
     t1 = ctet + cal*stet*sfi/sal;
-    assert(sal != 0);
     t2 = stet/sal;                                                       
     q1 = xq*(b1*t1 + b2*t2*cfi)/beta;
@@ -6220 +6142 @@
   hel = 0.0;                                                       
   w1 = q1*q1 + q2*q2 + q3*q3;
-  assert((w1 + x2*x2) >= 0);
   wq = std::sqrt(w1 + x2*x2);
   p1 = -q1;
   p2 = -q2;
   p3 = -q3;
-  assert((w1+x1*x1) >= 0);
   wp = std::sqrt(w1+x1*x1);
   loren(&q1, &q2, &q3, &b1, &b2, &b3, &wq);
@@ -6251 +6171 @@
 {
   // time 
-  G4double t[10];
+  G4double t[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
 
   t[0] = bl1->p1[i]/bl1->eps[i] - bl1->p1[j]/bl1->eps[j]; // t(1)->t[0] 
@@ -6275 +6195 @@
 void G4Incl::newt(G4int l1, G4int l2)
 {
-  G4int ig, id, kg, kd;
-  G4int iy, ix;
-  G4double E;
+  G4int ig = 0, id = 0, kg = 0, kd = 0;
+  G4int iy = 0, ix = 0;
+  G4double E = 0.0;
 
   G4int ia = bl3->ia1 + bl3->ia2;
@@ -6349 +6269 @@
     bl6->isa=bl1->ind2[ix];
   newt61:
-    // assert(isnan(E) == false);
-    // assert(isnan(iy) == false);
-    // assert(isnan(bl1->ind2[ig]) == false);
-    // assert(isnan(bl1->ind2[id]) == false);
     if ((31.*(bl3->rab2)) > totalCrossSection(E,iy,bl1->ind2[ig]+bl1->ind2[id])) {
       goto newt50;
@@ -6387 +6303 @@
     bl6->isa=bl1->ind2[ix];
   newt63: 
-    // assert(isnan(E) == false);
-    // assert(isnan(iy) == false);
-    // assert(isnan(bl1->ind2[kg]) == false);
-    // assert(isnan(bl1->ind2[kd]) == false);
     if ((31.*(bl3->rab2)) > totalCrossSection(E,iy,bl1->ind2[kg]+bl1->ind2[kd])) {
       continue;
@@ -6404 +6316 @@
 void G4Incl::new1(G4int l1)
 {
-  G4int ia, iy, ix;
-  G4double E;
+  G4int ia = 0, iy = 0, ix = 0;
+  G4double E = 0.0;
 
   ia=bl3->ia1+bl3->ia2;
@@ -6445 +6357 @@
     bl6->isa=bl1->ind2[ix];
   new161: 
-    // assert(isnan(E) == false);
-    // assert(isnan(iy) == false);
-    // assert(isnan(bl1->ind2[i]) == false);
-    // assert(isnan(bl1->ind2[l1]) == false);
     if ((31.*(bl3->rab2)) > totalCrossSection(E ,iy,bl1->ind2[i]+bl1->ind2[l1])) {
       continue;
@@ -6481 +6389 @@
     bl6->isa=bl1->ind2[ix];
   new171: 
-    // assert(isnan(E) == false);
-    // assert(isnan(iy) == false);
-    // assert(isnan(bl1->ind2[i]) == false);
-    // assert(isnan(bl1->ind2[l1]) == false);
     if ((31.0*(bl3->rab2)) > totalCrossSection(E,iy,bl1->ind2[i]+bl1->ind2[l1])) {
       continue;
@@ -6504 +6408 @@
                   G4double q4, G4int npion, G4int l1)
 {
-  G4double t[10];
+  G4double t[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
 
   G4int ia = bl3->ia1 + bl3->ia2;
@@ -6518 +6422 @@
     }
 
-    assert(q4 != 0);
-    assert(bl1->eps[i] != 0);
     t[0] = bl1->p1[i]/bl1->eps[i] - q1/q4;
     t[1] = bl1->p2[i]/bl1->eps[i] - q2/q4;
@@ -6539 +6441 @@
     G4double xx2 = t[3]*t[3] + t[4]*t[4] + t[5]*t[5] + (bl1->ta)*t[6];
     G4double E = std::sqrt(std::pow((bl1->eps[i]+q4),2) - std::pow((bl1->p1[i]+q1),2) - std::pow((bl1->p2[i]+q2),2) - std::pow((bl1->p3[i]+q3),2));
-    // assert(isnan(E) == false);
     if ((31.0*xx2) > pionNucleonCrossSection(E)) {
       continue;
@@ -6553 +6454 @@
                   G4double q4, G4int npion, G4int l1)
 {
-  G4double t[10];
-  G4double E, xx2;
-  G4int ia;
+  G4double t[10] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+  G4double E = 0.0, xx2 = 0.0;
+  G4int ia = 0;
 
   if(bl5->nesc[l1] > 0) {
@@ -6561 +6462 @@
   }
 
-  assert(q4 != 0);
-  assert(bl1->eps[l1] != 0);
   t[0] = bl1->p1[l1]/bl1->eps[l1] - q1/q4;
   t[1] = bl1->p2[l1]/bl1->eps[l1] - q2/q4;
@@ -6585 +6484 @@
   xx2 = t[3]*t[3] + t[4]*t[4] + t[5]*t[5] + (bl1->ta)*t[6];
   E = std::sqrt(std::pow((bl1->eps[l1]+q4),2) - std::pow((bl1->p1[l1]+q1),2) - std::pow((bl1->p2[l1]+q2),2) - std::pow((bl1->p3[l1]+q3),2));
-  // assert(isnan(E) == false);
   if ((31.0*xx2) > pionNucleonCrossSection(E)) {
     return;
@@ -6633 +6531 @@
   //   3760           common/ws/r0,adif,rmaxws,drws,nosurf,xfoisa,npaulstr,bmax
 
-  G4double pmod, pr2;
-  G4double xr2, rdeq, dx2, dp2;
-  G4double rs, vol;
-  G4int nl;
-  G4int ia;
+  G4double pmod = 0.0, pr2 = 0.0;
+  G4double xr2 = 0.0, rdeq = 0.0, dx2 = 0.0, dp2 = 0.0;
+  G4double rs = 0.0, vol = 0.0;
+  G4int nl = 0;
+  G4int ia = 0;
 
   if (ws->npaulstr == 2) {
@@ -6648 +6546 @@
     if (pmod < 270.0) {
       return 1.0;
-    }
-    else {
+    } else {
       return 0.0;
     }
-  }
-  else {
+  } else {
     // Statistic Pauli blocking
     xr2 = xr*xr;
@@ -6659 +6555 @@
     vol = std::pow((40.0*3.1415926/3.0),2) * (std::pow((xr*pr)/(2.0*3.1415926*197.13),3));
     rs = std::sqrt(bl3->x1[l]*bl3->x1[l] + bl3->x2[l]*bl3->x2[l] + bl3->x3[l]*bl3->x3[l]);
-    // assert(isnan(rs) == false);
     if (ws->nosurf <= 0) {
       // modifs a.b.: r2 -> rmaxws pour la densite en w.s.
       rdeq = ws->rmaxws;
-    }
-    else {
+    } else {
       rdeq = ws->r0;
     }
@@ -6683 +6577 @@
           if(((nl - 1)/vol/2.0) > 1.0) {
             return 1.0;
-          }
-          else {
+          } else {
             return ((nl - 1)/vol/2.0);
           }
@@ -6690 +6583 @@
         nl = nl + 1;
       }
-    }
-    else {
+    } else {
       return 0.0;
     }
@@ -6709 +6601 @@
   G4double plab = E*std::sqrt(E*E-3.52e6)/1876.6;
   G4double p1 = 0.001*plab;
-  G4double alp;
+  G4double alp = 0.0;
 
   if(plab > 2000.0) {
@@ -6756 +6648 @@
   // m=0,1,2 for nucleon-nucleon,nucleon-delta,delta,delta
 
-  G4double stotResult;
+  G4double stotResult = 0.0;
   G4double sine = 0.0;
 
   if((m-1) < 0) {
-    // assert(isnan(E) == false);
-    // assert(isnan(i) == false);
     sine = deltaProductionCrossSection(E,int(i));
   }
@@ -6774 +6664 @@
 
   stotResult = sine + lowEnergy(E,m,i);
-  // assert(isnan(stotResult) == false);
   return stotResult;
 }
@@ -6781 +6670 @@
 {
   G4double E = Ein;
-  G4double s;
-  G4double x, y;
-  G4double srecResult;
+  G4double s = 0.0;
+  G4double x = 0.0, y = 0.0;
+  G4double srecResult = 0.0;
 
   if (i*i == 16) {
@@ -6791 +6680 @@
   if(E <= (938.3 + d)) {
     return 0.0;
-  }
-  else {
+  } else {
     if(E < (938.3 + d + 2.0)) {
       E = 938.3 + d + 2.0;
@@ -6799 +6687 @@
     x = (s - 3.523e6)/(s - std::pow((938.3 + d),2));
     y = s/(s - std::pow((d - 938.3),2));
-    // assert(isnan(E) == false);
-    // assert(isnan(i) == false);
     srecResult = 0.5*x*y*deltaProductionCrossSection(E, i);
     srecResult = srecResult*(32.0 + i*i*(isa*isa - 5))/64.0;
@@ -6877 +6763 @@
 //   }
   double scali=1.0;
-  double plab;
-  double p1;
-  double sproResult;
-
-  double EE=E-(bl8->rathr);
-
-  // assert(isnan(EE) == false);
+  double plab = 0.0;
+  double p1 = 0.0;
+  double sproResult = 0.0;
+
+  double EE = E -(bl8->rathr);
+
   if(EE*EE-3.53e6 < 0) {
     goto spro22;
   }
-  plab=EE*std::sqrt(EE*EE-3.52e6)/1876.6;
-  // assert(isnan(plab) == false);
-  p1=0.001*plab;
+  plab = EE * std::sqrt(EE*EE - 3.52e6)/1876.6;
+  p1 = 0.001*plab;
   if (plab > 800.) {
     goto spro1;
   }
   spro22:
-  sproResult=0.0;
+  sproResult = 0.0;
   return sproResult;
 
@@ -6911 +6795 @@
   sproResult=(24.2+8.9*p1-31.1/std::sqrt(p1))*scali;
   return sproResult;
-  spro3: sproResult=(33.0 + 196.0*std::sqrt(std::pow(std::fabs(p1-0.95),5))-31.1/std::sqrt(p1))*scali;
+  spro3: sproResult = (33.0 + 196.0*std::sqrt(std::pow(std::fabs(p1 - 0.95),5)) - 31.1/std::sqrt(p1))*scali;
   return sproResult;
   spro10: if (plab < 2000.) {
     goto spro11;
   }
-  sproResult=(41.+(60.*p1-54.)*std::exp(-1.2*p1)-77./(p1+1.5))*scali;
+  sproResult = (41.0 + (60.0*p1 - 54.0)*std::exp(-1.2*p1) - 77.0/(p1 + 1.5))*scali;
   return sproResult;
   spro11: if (plab < 1500.) {
     goto spro12;
   }
-  sproResult=41.+60.*(p1-0.9)*std::exp(-1.2*p1)-1250./(p1+50.)+4.*std::pow((p1-1.3),2);
+  sproResult = 41.0 + 60.0*(p1 - 0.9)*std::exp(-1.2*p1) - 1250.0/(p1 + 50.0)+4.0*std::pow((p1-1.3),2);
+  sproResult = sproResult*scali;
+  return sproResult;
+ spro12: sproResult=23.5 + 24.6/(1.0 + std::exp(-10.0*p1 + 12.0)) - 1250.0/(p1 + 50.0) + 4.0*std::pow((p1 - 1.3),2);
   sproResult=sproResult*scali;
   return sproResult;
- spro12: sproResult=23.5+24.6/(1.+std::exp(-10.*p1+12.))-1250./(p1+50.)+4.*std::pow((p1-1.3),2);
-  sproResult=sproResult*scali;
-  return sproResult;
-
 }
 
@@ -6938 +6821 @@
   G4double y = x*x;
   G4double q2 = (y-std::pow(1076.0,2))*(y-std::pow(800.0,2))/y/4.0;
-  assert(q2 >= 0);
-  G4double q3, f3;
-  G4double spn;
+  G4double q3 = 0.0, f3 = 0.0;
+  G4double spn = 0.0;
 
   if(q2 <= 0) {
     return 0.0;
-  }
-  else {
+  } else {
     q3 = std::pow((std::sqrt(q2)),3);
     f3 = q3/(q3+std::pow(180.0,3));
@@ -6961 +6842 @@
   // of the nucleus and r is the target radius
 
-  G4double x;
+  G4double x = 0.0;
   G4double barr = 0.0;
 
   // We need enough energy to escape from the potential well.
   if (E > v0) {
-    x = std::sqrt(E*(E-v0));
-    // assert(isnan(x) == false);
-    barr = 4.*x/(E+E-v0+x+x);
-    // assert(isnan(barr) == false);
+    x = std::sqrt(E*(E - v0));
+    barr = 4.0*x/(E + E - v0 + x + x);
     if (iz > 0) {
       G4double b = izn*1.44/r;
-      G4double px = std::sqrt((E-v0)/b);
-      // assert(isnan(px) == false);
+      G4double px = std::sqrt((E - v0)/b);
       
       if (px < 1.0) {
-        G4double g = izn/137.03*std::sqrt(2.*938.3/(E-v0))*(std::acos(px)-px*std::sqrt(1.-px*px));
-        // assert(isnan(g) == false);
+        G4double g = izn/137.03*std::sqrt(2.0*938.3/(E - v0))*(std::acos(px) - px*std::sqrt(1.0 - px*px));
         if (g > 35.){
           barr=0.0;
-        }
-        else {
+        } else {
           barr = barr*std::exp(-2.0*g);
         }
         return barr;
-      }
-      else {
+      } else {
         return barr;
       }
-    }
-    else {
+    } else {
       return barr;
     }
-  }
-  else {
+  } else {
     return barr;
   }
@@ -7003 +6876 @@
   const G4double  pf = 270.339 , pf2 = 73083.4;
 
-  G4double ref;
-  G4double t1, t3, t4, t5;
+  G4double ref = 0.0;
+  G4double t1 = 0.0, t3 = 0.0, t4 = 0.0, t5 = 0.0;
   
-  G4double t2 = p1*p1+p2*p2+p3*p3;
-  assert(t2 >= 0);
+  G4double t2 = p1*p1 + p2*p2 + p3*p3;
   G4double p = std::sqrt(t2);
   G4double r = r2;
-  G4double xv;
-  G4double s;
+  G4double xv = 0.0;
+  G4double s = 0.0;
   
   if (ws->nosurf <= 0) {
     xv = p/pf;
-    // assert(isnan(xv) == false);
     r = interpolateFunction(xv);
-    // assert(isnan(r) == false);
     r = r*r;
     if (t2 > pf2) {
@@ -7039 +6909 @@
   if (t5 > 0) {
     ref = (-t3 + std::sqrt(t5))*E;
-    // assert(isnan(ref) == false);
     return ref;
   }
@@ -7078 +6947 @@
   }
 
-  bl3->ia2 = int(calincl->f[0]); // f(1) -> f[0]
+  bl3->ia2 = int(std::floor(calincl->f[0] + 0.1)); // f(1) -> f[0]
   sep = 6.8309;
 
@@ -7105 +6974 @@
       bmaxt = ws->bmax;
       proba_trans = coulombTransm(calincl->f[2],bl3->ia1,iz1,calincl->f[0],calincl->f[1]);
-      // assert(isnan(proba_trans) == false);
-      
       proba = forceAbs(1,calincl->f[0],calincl->f[1],calincl->f[2],bmaxt,proba_trans);
-      // assert(isnan(proba) == false);
       
       standardRandom(&alea,&(hazard->igraine[4]));
@@ -7119 +6985 @@
       
       del = std::sqrt(std::pow(((calincl->f[0] + 1.0)*fmpinc + calincl->f[2]),2) - pbeam2);
-      // assert(isnan(del) == false);
       
       (*erecrem) = pbeam2/((calincl->f[0] + 1.0)*fmpinc+calincl->f[2] + del);
@@ -7139 +7004 @@
       bmaxt = ws->bmax;
       proba_trans = coulombTransm(calincl->f[2],bl3->ia1,iz1,calincl->f[0],calincl->f[1]);
-      // assert(isnan(proba_trans) == false);
-      
       proba = forceAbs(1,calincl->f[0],calincl->f[1],calincl->f[2],bmaxt,proba_trans);
-      // assert(isnan(proba) == false);
       
       standardRandom(&alea,&(hazard->igraine[4]));
@@ -7153 +7015 @@
       
       del = std::sqrt(std::pow(((calincl->f[0]+1.)*fmpinc+calincl->f[2]),2)-pbeam2);
-      // assert(isnan(del) == false);
       
       (*erecrem) = pbeam2/((calincl->f[0] + 1.0)*fmpinc + calincl->f[2] + del);
@@ -7190 +7051 @@
 
   sig_exp = xabs2(zp, ap, zt, at, ep);
-  // assert(isnan(sig_exp) == false);
-  
   sig_incl = crossSection(int(iprojo), ep, at);
-  // assert(isnan(sig_incl) == false);
   
   proba = (sig_exp-pt*sig_incl)/(pt*(sig_g - sig_incl));
@@ -7268 +7126 @@
   gcm = (ap*(dp1 + E/938.0) + at)/(std::pow((std::pow(ap,2) + std::pow(at,2) + dp2*ap*(E + 938.0)*at/938.e0),dph));
   bcm = std::sqrt(dp1-dp1/(std::pow(gcm,2)));
-  // assert(isnan(bcm) == false);
   
   plab = ap*std::sqrt(dp2*938.0*E + E*E);                                    
-  // assert(isnan(plab) == false);
   ecmp = gcm*(E+938.0)*ap - bcm*gcm*plab - ap*938.0;                     
   ecmt = gcm*938.0*at - at*938.0;
@@ -7281 +7137 @@
   rm = (197.32/137.01)*zp*zt/ecm;                                      
   bigr = rp + rt + 1.2*(std::pow(ap,dpth) + std::pow(at,dpth))/(std::pow(ecm,dpth));                   
-  // assert(isnan(bigr) == false);
 
   bigb = 1.44*zp*zt/bigr;                                              
-  // assert(isnan(bigb) == false);
   
   if ((nint(zp) == 1) && (nint(at) > 56)) {
@@ -7420 +7274 @@
   }
 
-//  if(isnan(sig)) {
-//    sig = 0.0;
-//  }
-
   return sig;
 }
@@ -7431 +7281 @@
   (*seed) = (*seed); // Avoid warning during compilation.
   // Use Geant4 G4UniformRand
-  (*rndm) = G4UniformRand();
+  //  (*rndm) = G4UniformRand();
+  (*rndm) = randomGenerator->getRandom();
 }
 
@@ -7516 +7367 @@
                                  {-0.01087e+0,2.6494e+0,-2.5173e+0}}; 
 
-  G4double apow[3], epow[5];
-  G4int ii, jj;
+  G4double apow[3] = {0.0, 0.0, 0.0}, epow[5] = {0.0, 0.0, 0.0};
+  G4int ii = 0, jj = 0;
 
   if(A >= 27.0) {
     ii = 3;        
     jj = 4;
-  }
-  else {
+  } else {
     ii = 3;
     jj = 5;
@@ -7545 +7395 @@
         }
       }
-    }
-    else {
+    } else {
       for(G4int i = 0; i < ii; i++ ) {
         for(G4int j = 0; j < jj; j++) {
@@ -7553 +7402 @@
       }
     }
-  }
-  else {
+  } else {
     for(G4int i = 0; i < ii; i++) {
       for(G4int j = 0; j < jj; j++) {
@@ -7567 +7415 @@
 G4double G4Incl::coulombTransm(G4double E, G4double fm1, G4double z1, G4double fm2, G4double z2)
 {
-  G4double eta,rho;
+  G4double eta = 0.0, rho = 0.0;
   const G4double c2 = 0.00516;
   const G4double c3 = 0.007165;
   const G4double uma = 938.0;
-  G4double ml;
+  G4double ml = 0.0;
   
   G4double ecm = E*fm2/(fm1+fm2);
@@ -7603 +7451 @@
   static G4double psi = rho*y;                                                         
 
-  static G4int i0, j0;
+  static G4int i0 = 0, j0 = 0;
   
-  static G4double prob;
-  static G4double dumm, x, cx;
-  static G4double t, t1, t2, t3;
-  static G4double f, g;
-  static G4double temp, temp1, temp2;
-  static G4double xk, delp0, delp1, delx0, delx1;
+  static G4double prob = 0.0;
+  static G4double dumm = 0.0, x = 0.0, cx = 0.0;
+  static G4double t = 0.0, t1 = 0.0, t2 = 0.0, t3 = 0.0;
+  static G4double f = 0.0, g = 0.0;
+  static G4double temp = 0.0, temp1 = 0.0, temp2 = 0.0;
+  static G4double xk = 0.0, delp0 = 0.0, delp1 = 0.0, delx0 = 0.0, delx1 = 0.0;
 
   if (rho > y) {                                                 
     if (psi > dp4 && psi < 50.0) {                           
       prob = clmb2(rho,eta,&dumm);                                       
-    }
-    else {
+    } else {
       x = std::exp(std::log(eta)/6.0);
       prob = std::sqrt(dp1 - y*x/(c0 + c1 * std::pow(x,3) + rho * x));
     } 
     (*ml) = 0;
-  }                                                            
-  else {
+  } else {
     x = rho/y;                                                          
     if (psi <= psi0[0]) {
@@ -7631 +7477 @@
       if (eta > dp1) {                                          
         t3 = std::log(psi)+dp2*gamma - dp1 + dp1/(12.e0*std::pow(eta,2))+dp1/(12.e1*std::pow(eta,4));
-      }                                                        
-      else {                                                          
+      } else {                                                          
         t3 = std::log(dp2*rho) + gamma - dp1/(dp1 + std::pow(eta,2)) + s3*std::pow(eta,2) + s4*std::pow(eta,4);   
       }
@@ -7639 +7484 @@
       prob = cx/(std::pow(g,2)+std::pow(f,2));                                           
       (*ml) = 3;
-    }                                                          
-    else if (psi <= psi0[ln0-1]) {
+    } else if (psi <= psi0[ln0-1]) {
       if (x <= x0[0]) {
         temp = std::log(psi/psi0[0]);
@@ -7653 +7497 @@
         prob = cx/std::pow(prob,2);                                             
         (*ml) = 1;
-      }                                                        
-      else {                                                          
+      } else {                                                          
         temp1 = std::log(x/x0[0]);
         i0 = min(max(1 + int(temp1/delx0),1),lt0-1);                                     
@@ -7666 +7509 @@
         (*ml)=2;                                                        
       }
-    }
-    else if (psi <= psi1[ln1-1]) {
+    } else if (psi <= psi1[ln1-1]) {
       if (x <= x1[0]) {
         temp = std::log(psi/psi1[0]);
@@ -7678 +7520 @@
         prob = cx/std::pow(prob,2);                                             
         (*ml) = 1;
-      }                                                        
-      else {                                                         
+      } else {                                                         
         temp1 = std::log(x/x1[0]);
         i0 = min(max(1+int(temp1/delx1),1),lt1-1);                                    
@@ -7691 +7532 @@
         (*ml)=2;                                                        
       }
-    }
-    else {
+    } else {
       prob = clmb2(rho,eta,&dumm);                                       
       (*ml) = 4;                                              
@@ -7722 +7562 @@
   const G4double x1 = 0.01;
   const G4double xi = 100;                                        
-  static G4double x,temp,prob;
-  static G4int i;
+  static G4double x = 0.0, temp = 0.0, prob = 0.0;
+  static G4int i = 0;
 
   x = dp1/(dp1 + std::sqrt(dph*rho*eta));                                     
   if (x < x1) {                                                 
     temp = t0[2] * std::pow((x/x1),dpth);
-  }                                       
-  else {
+  } else {
     i = int(std::floor(xi*x));                                                          
     i = i + 1;
@@ -7750 +7589 @@
   if(a < b) {
     return a;
-  }
-  else {
+  } else {
     return b;
   }
@@ -7760 +7598 @@
   if(a < b) {
     return a;
-  }
-  else {
+  } else {
     return b;
   }
@@ -7770 +7607 @@
   if(a > b) {
     return a;
-  }
-  else {
+  } else {
     return b;
   }
@@ -7780 +7616 @@
   if(a > b) {
     return a;
-  }
-  else {
+  } else {
     return b;
   }
@@ -7797 +7632 @@
     if(fractpart < 0.5) {
       return int(std::floor(number));
-    }
-    else {
+    } else {
       return int(std::ceil(number));
     }
@@ -7805 +7639 @@
     if(fractpart < -0.5) {
       return int(std::floor(number));
-    }
-    else {
+    } else {
       return int(std::ceil(number));
     }
@@ -7845 +7678 @@
 G4double G4Incl::pcm(G4double E, G4double A, G4double C)
 {
-//   assert(((std::pow(E,2)-std::pow((A+C),2))*(std::pow(E,2)-std::pow((A-C),2))) >= 0);
-//   assert(E != 0);
   return (0.5*std::sqrt((std::pow(E,2)-std::pow((A+C),2))*(std::pow(E,2)-std::pow((A-C),2)))/E);
 }

trunk/source/processes/hadronic/models/incl/src/G4InclAblaDataFile.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4InclAblaDataFile.cc,v 1.5 2007/12/03 19:36:06 miheikki Exp $ 
+// $Id: G4InclAblaDataFile.cc,v 1.7 2008/06/25 17:20:04 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -87 +87 @@
   
   G4double flalpha, frldm, vgsld, pace2;
-  const G4int rows = 98;
-  const G4int cols = 153;
+  const G4int rows = 99;
+  const G4int cols = 154;
   const G4int massnumbers = 263;
-  for(int i = 0; i < rows; i++) {
-    for(int j = 0; j < cols; j++) {
+  for(int i = 0; i < cols; i++) {
+    for(int j = 0; j < rows; j++) {
       setAlpha(j, i, 0.0);
       setEcnz( j, i, 0.0);
@@ -98 +98 @@
   }
   
-  for(int i = 0; i < rows; i++) {
-    for(int j = 0; j < cols; j++) {
+  for(int i = 0; i < cols; i++) {
+    for(int j = 0; j < rows; j++) {
       flalphain >> flalpha;
       frldmin >> frldm;
       vgsldin >> vgsld;      
-      setAlpha(j, i, flalpha);
-      setEcnz( j, i, frldm);
-      setVgsld(j, i, vgsld);      
+      setAlpha(i, j, flalpha);
+      setEcnz( i, j, frldm);
+      setVgsld(i, j, vgsld);      
     }
   }