Index: trunk/source/processes/hadronic/models/incl/src/G4Abla.cc
===================================================================
--- trunk/source/processes/hadronic/models/incl/src/G4Abla.cc	(revision 1340)
+++ trunk/source/processes/hadronic/models/incl/src/G4Abla.cc	(revision 1347)
@@ -24,5 +24,5 @@
 // ********************************************************************
 //
-// $Id: G4Abla.cc,v 1.22 2010/10/26 02:47:59 kaitanie Exp $ 
+// $Id: G4Abla.cc,v 1.27 2010/11/17 20:19:09 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -40,31 +40,4 @@
 #include "G4AblaFissionSimfis18.hh"
 #include "G4AblaFission.hh"
-
-G4Abla::G4Abla()
-{
-  verboseLevel = 0;
-  ilast = 0;
-}
-
-G4Abla::G4Abla(G4Hazard *hazard, G4Volant *volant)
-{
-  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();
-  eenuc = new G4Eenuc();
-  ec2sub = new G4Ec2sub();
-  ecld = new G4Ecld();
-  fb = new G4Fb();
-  fiss = new G4Fiss();
-  opt = new G4Opt();
-}
 
 G4Abla::G4Abla(G4Hazard *aHazard, G4Volant *aVolant, G4VarNtp *aVarntp)
@@ -109,4 +82,5 @@
 G4Abla::~G4Abla()
 {
+  delete fissionModel;
   delete randomGenerator;
   delete pace;
@@ -483,6 +457,6 @@
       // 	       G4double *pleva_par, G4double *pxeva_par, G4double *pyeva_par,
       // 	       G4double *ff_par, G4int *inttype_par, G4int *inum_par);
-      G4double zf1, af1, malpha1, ffpleva1, ffpxeva1, ffpyeva1;
-      G4int ff1, ftype1;
+      G4double zf1 = 0.0, af1 = 0.0, malpha1 = 0.0, ffpleva1 = 0.0, ffpxeva1 = 0.0, ffpyeva1 = 0.0;
+      G4int ff1 = 0, ftype1 = 0;
       evapora(zff1, aff1, epf1_out, 0.0, &zf1, &af1, &malpha1, &ffpleva1,
 	      &ffpxeva1, &ffpyeva1, &ff1, &ftype1, &inum);
@@ -560,6 +534,6 @@
       // 	       G4double *pleva_par, G4double *pxeva_par, G4double *pyeva_par,
       // 	       G4double *ff_par, G4int *inttype_par, G4int *inum_par);
-      G4double zf2, af2, malpha2, ffpleva2, ffpxeva2, ffpyeva2;
-      G4int ff2, ftype2;
+      G4double zf2 = 0.0, af2 = 0.0, malpha2 = 0.0, ffpleva2 = 0.0, ffpxeva2 = 0.0, ffpyeva2 = 0.0;
+      G4int ff2 = 0, ftype2 = 0;
       evapora(zff2,aff2,epf2_out,0.0,&zf2,&af2,&malpha2,&ffpleva2,
 	      &ffpxeva2,&ffpyeva2,&ff2,&ftype2,&inum);
@@ -1608,5 +1582,5 @@
 		    G4double *sbp_par, G4double *sba_par, G4double *ecn_par, 
 		    G4double *ecp_par,G4double *eca_par, G4double *bp_par,
-		    G4double *ba_par, G4int inttype, G4int inum, G4int itest)
+		    G4double *ba_par, G4int, G4int inum, G4int itest)
 {
   G4int dummy0 = 0;
@@ -1758,5 +1732,4 @@
 
   imaxwell = 1;
-  inttype = 0;
   
   // limiting of excitation energy where fission occurs                    
@@ -2842,7 +2815,7 @@
   i = idint(y/(2.0e-02)) + 1;
     
-  if(i >= bsbkSize) {
+  if((i + 1) >= bsbkSize) {
     if(verboseLevel > 2) {
-      G4cout <<"G4Abla error: index i = " << i << " is greater than array size permits." << G4endl;
+      G4cout <<"G4Abla error: index " << i + 1 << " is greater than array size permits." << G4endl;
     }
     bipolResult = 0.0;
@@ -3704,7 +3677,6 @@
 // own class
 
-void G4Abla::standardRandom(G4double *rndm, G4long *seed)
-{
-  (*seed) = (*seed); // Avoid warning during compilation.
+void G4Abla::standardRandom(G4double *rndm, G4long*)
+{
   // Use Geant4 G4UniformRand
   (*rndm) = randomGenerator->getRandom();
Index: trunk/source/processes/hadronic/models/incl/src/G4AblaFission.cc
===================================================================
--- trunk/source/processes/hadronic/models/incl/src/G4AblaFission.cc	(revision 1340)
+++ trunk/source/processes/hadronic/models/incl/src/G4AblaFission.cc	(revision 1347)
@@ -24,5 +24,5 @@
 // ********************************************************************
 //
-// $Id: G4AblaFission.cc,v 1.3 2008/11/06 08:42:00 gcosmo Exp $
+// $Id: G4AblaFission.cc,v 1.5 2010/11/17 20:19:09 kaitanie Exp $
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -36,4 +36,6 @@
 G4AblaFission::G4AblaFission()
 {
+  hazard = 0;
+  randomGenerator = 0;
 }
 
@@ -1086,7 +1088,6 @@
 }
 
-void G4AblaFission::standardRandom(G4double *rndm, G4long *seed)
-{
-  (*seed) = (*seed); // Avoid warning during compilation.
+void G4AblaFission::standardRandom(G4double *rndm, G4long*)
+{
   // Use Geant4 G4UniformRand
   (*rndm) = randomGenerator->getRandom();
Index: trunk/source/processes/hadronic/models/incl/src/G4AblaFissionSimfis18.cc
===================================================================
--- trunk/source/processes/hadronic/models/incl/src/G4AblaFissionSimfis18.cc	(revision 1340)
+++ trunk/source/processes/hadronic/models/incl/src/G4AblaFissionSimfis18.cc	(revision 1347)
@@ -24,5 +24,5 @@
 // ********************************************************************
 //
-// $Id: G4AblaFissionSimfis18.cc,v 1.3 2008/11/06 08:42:00 gcosmo Exp $
+// $Id: G4AblaFissionSimfis18.cc,v 1.5 2010/11/17 20:19:09 kaitanie Exp $
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -36,4 +36,6 @@
 G4AblaFissionSimfis18::G4AblaFissionSimfis18()
 {
+  hazard = 0;
+  randomGenerator = 0;
 }
 
@@ -1412,7 +1414,6 @@
 }
 
-void G4AblaFissionSimfis18::standardRandom(G4double *rndm, G4long *seed)
-{
-  (*seed) = (*seed); // Avoid warning during compilation.
+void G4AblaFissionSimfis18::standardRandom(G4double *rndm, G4long*)
+{
   // Use Geant4 G4UniformRand
   (*rndm) = randomGenerator->getRandom();
Index: trunk/source/processes/hadronic/models/incl/src/G4Incl.cc
===================================================================
--- trunk/source/processes/hadronic/models/incl/src/G4Incl.cc	(revision 1340)
+++ trunk/source/processes/hadronic/models/incl/src/G4Incl.cc	(revision 1347)
@@ -24,5 +24,5 @@
 // ********************************************************************
 //
-// $Id: G4Incl.cc,v 1.33 2010/10/29 06:48:43 gunter Exp $ 
+// $Id: G4Incl.cc,v 1.37 2010/12/15 07:41:31 gunter Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -57,4 +57,5 @@
   randomGenerator = new G4InclGeant4Random();
   //  randomGenerator = new G4Ranecu();
+  clearState();
 }
 
@@ -85,4 +86,5 @@
   //  randomGenerator = new G4Ranecu();
   randomGenerator = new G4InclGeant4Random();
+  clearState();
 }
 
@@ -149,8 +151,14 @@
 
   theLogger = 0;
+  clearState();
 }
 
 G4Incl::~G4Incl()
 {
+  delete be;
+  delete ps;
+  delete fermi;
+  delete qvp;
+
   delete randomGenerator;
   delete light_gaus_nuc;
@@ -385,4 +393,5 @@
 void G4Incl::processEventIncl(G4InclInput *input)
 {
+  //G4cout <<"Starting event " << eventnumber << G4endl;
   if(input == 0) {
     G4cerr <<"G4Incl fatal error: NULL pointer passed as input!" << G4endl;
@@ -393,14 +402,122 @@
   const G4double uma = 931.4942;
   const G4double melec = 0.511;
-
-  G4double   pcorem = 0.0;
-  G4double   pxrem = 0.0;
-  G4double   pyrem = 0.0;
-  G4double   pzrem = 0.0;
+  const G4double fmp = 938.2796;
+
+  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->bulletType() == 3) { // pi+ 
+
+  if(calincl->bulletType() == -12) {
+    be->ia_be = std::abs(calincl->bulletType());
+    be->iz_be = 6;
+  } else if(calincl->bulletType() == -666) {
+    be->iz_be = calincl->extendedProjectileZ();
+    be->ia_be = calincl->extendedProjectileA();
+  }
+
+  if(calincl->isExtendedProjectile() == false && calincl->bulletType() < -max_a_proj) {
+  //  if(calincl->bulletType() < -max_a_proj) {
+    G4cout <<"max a of composite projectile is: " << max_a_proj << G4endl;
+    exit(0);
+  }
+  if(calincl->bulletType() < 0) {
+    //    calincl->bulletType() = std::floor(calincl->bulletType() + 0.1); WTF???
+    be->pms_be=100.;
+    G4int i_tabled=0;
+    if(be->iz_be == 3 && be->ia_be == 6) {
+      be->rms_be=2.56;
+      be->bind_be=32.0;
+      i_tabled=1;
+    } else if(be->iz_be == 3 && be->ia_be == 7) { // TODO: Check the values!
+      be->rms_be=2.56;
+      be->bind_be=32.0;
+      i_tabled=1;
+    } else if(be->iz_be == 3 && be->ia_be == 8) {
+      be->rms_be=2.40;
+      be->bind_be=39.25;
+      i_tabled=1;
+    } else if(be->iz_be == 4 && be->ia_be == 7) {
+      be->rms_be=2.51;
+      be->bind_be=58.17;
+      i_tabled=1;
+    } else if(be->iz_be == 4 && be->ia_be == 9) {
+      be->rms_be=2.51;
+      be->bind_be=58.17;
+      i_tabled=1;
+    } else if(be->iz_be == 4 && be->ia_be == 10) {
+      be->rms_be=2.45;
+      be->bind_be=64.75;
+      i_tabled=1;
+    } else if(be->iz_be == 5 && be->ia_be == 10) {
+      be->rms_be=2.45;
+      be->bind_be=64.75;
+      i_tabled=1;
+    } else if(be->iz_be == 5 && be->ia_be == 11) {
+      be->rms_be=2.40;
+      be->bind_be=76.21;
+      i_tabled=1;
+    } else if(be->iz_be == 6 && be->ia_be == 9) { // TODO: Check the values!
+      be->rms_be=2.44;
+      be->bind_be=92.17;
+      i_tabled=1;
+    } else if(be->iz_be == 6 && be->ia_be == 10) { // TODO: Check the values!
+      be->rms_be=2.44;
+      be->bind_be=92.17;
+      i_tabled=1;
+    } else if(be->iz_be == 6 && be->ia_be == 11) { // TODO: Check the values!
+      be->rms_be=2.44;
+      be->bind_be=92.17;
+      i_tabled=1;
+    } else if(be->iz_be == 6 && calincl->bulletType() == -12) { // Special Carbon case
+      G4cout <<"Carbon 12 (special) selected." << G4endl;
+      be->rms_be=2.44;
+      be->bind_be=92.17;
+      i_tabled=1;
+    } else if(be->iz_be == 6 && be->ia_be == 12) {
+      be->rms_be=2.44;
+      be->bind_be=92.17;
+      i_tabled=1;
+    } else if(be->iz_be == 7 && be->ia_be == 16) {
+      be->rms_be=2.73;
+      be->bind_be=127.62;
+      i_tabled=1;
+    } else {
+      G4cout <<"Warning: No rms and binding for projectile ion A = " << be->ia_be << " Z = " << be->iz_be << G4endl;
+      be->rms_be=2.44;
+      be->bind_be=92.17;
+      G4cout <<"Warning: Using probably bad values rms = " << be->rms_be << " binding = " << be->bind_be << G4endl;
+      i_tabled=1;      
+    }
+      
+    if(i_tabled == 0) {
+      G4cout <<"This heavy ion (a,z)= " << be->ia_be << " " << be->iz_be << " is not defined as beam in INCL" << G4endl;
+      exit(0);
+    }
+      
+    //    G4cout <<"z projectile, rms_r, rms_p (gaussian model)" << be->iz_be << " " << be->rms_be << " " << be->pms_be << G4endl;
+    //    G4cout <<"binding energy (mev):" << be->bind_be << G4endl;
+    //    G4cout <<"fermi-breakup dresner below a=" << calincl->f[11] << G4endl;
+  }      
+  //  G4cout <<"Target Mass and Charge: " << calincl->targetA() << " " << calincl->targetZ() << G4endl;
+  //  calincl->f[10] = 0; // No clusters
+
+  if(calincl->bulletType() == -12) {  // C12 special case
+    mprojo=fmp*std::abs(calincl->bulletType()) - be->bind_be;
+    pbeam=std::sqrt(calincl->bulletE()*(calincl->bulletE()+2.*mprojo));
+    ap=std::abs(calincl->bulletType());
+    zp=be->iz_be;
+  } else if(calincl->bulletType() == -666) { // Generic extended projectile
+    mprojo=fmp*be->ia_be - be->bind_be;
+    pbeam=std::sqrt(calincl->bulletE()*(calincl->bulletE()+2.*mprojo));
+    ap=be->ia_be;
+    zp=be->iz_be;
+  }
+  // pi+
+  if(calincl->bulletType() == 3) { 
     mprojo = 139.56995;
     ap = 0.0;
@@ -408,5 +525,6 @@
   }
 
-  if(calincl->bulletType() == 4) { // pi0
+  // pi0
+  if(calincl->bulletType() == 4) {
     mprojo = 134.9764;
     ap = 0.0;
@@ -414,5 +532,6 @@
   }
 
-  if(calincl->bulletType() == 5) { // pi-
+  // pi-
+  if(calincl->bulletType() == 5) {
     mprojo = 139.56995;
     ap = 0.0;
@@ -420,30 +539,36 @@
   }
 
-  // Coulomb en entree seulement pour les particules ci-dessous.
-  if(calincl->bulletType() == 1) { // proton
+  // coulomb en entree seulement pour les particules ci-dessous
+
+  // proton
+  if(calincl->bulletType() == 1) {
     mprojo = 938.27231;
     ap = 1.0;
     zp = 1.0;
   }
-  
-  if(calincl->bulletType() == 2) { // neutron  
+
+  // neutron  
+  if(calincl->bulletType() == 2) {
     mprojo = 939.56563;
     ap = 1.0;
     zp = 0.0;
   }
-  
-  if(calincl->bulletType() == 6) { // deuteron
+
+  // deuteron
+  if(calincl->bulletType() == 6) {
     mprojo = 1875.61276;
     ap = 2.0;
     zp = 1.0;
   }
-  
-  if(calincl->bulletType() == 7) { // triton
+
+  // triton
+  if(calincl->bulletType() == 7) {
     mprojo = 2808.95;
     ap = 3.0;
     zp = 1.0;
   }
-  
-  if(calincl->bulletType() == 8) { // He3
+
+  // He3
+  if(calincl->bulletType() == 8) {
     mprojo = 2808.42;
     ap = 3.0;
@@ -451,5 +576,6 @@
   }
 
-  if(calincl->bulletType() == 9) { // Alpha
+  // Alpha
+  if(calincl->bulletType() == 9) {
     mprojo = 3727.42;
     ap = 4.0;
@@ -457,12 +583,9 @@
   }
 
-  if(calincl->bulletType() == -12) { // Carbon
-    mprojo = 12.0 * 938.27231;
+  // Carbon
+  if(calincl->bulletType() == -12) {
+    mprojo = 6.0*938.27231 + 6.0*939.56563;
     ap = 12.0;
     zp = 6.0;
-  } else if(calincl->bulletType() == -666) { // Generic extended ion projectile
-    mprojo = calincl->extendedProjectileA() * 938.27231;
-    ap = calincl->extendedProjectileA();
-    zp = calincl->extendedProjectileZ();
   }
 
@@ -482,13 +605,14 @@
   G4double bimpac = 0.0;
   G4int jrem = 0;
+  G4double xjrem = 0.0, yjrem = 0.0, zjrem = 0.0;
   G4double alrem = 0.0;
 
-  /**
-   * Coulomb barrier treatment.
-   */
+  // Coulomb barrier
+  
   G4double probaTrans = 0.0;
   G4double rndm = 0.0;
 
   if((calincl->bulletType() == 1) || (calincl->bulletType() >= 6)) {
+    //    probaTrans = coulombTransm(calincl->bulletE(),apro,zpro,calincl->targetA(),calincl->targetZ());
     probaTrans = coulombTransm(calincl->bulletE(),ap,zp,calincl->targetA(),calincl->targetZ());
     standardRandom(&rndm, &(hazard->ial));
@@ -499,33 +623,56 @@
   }
 
-  /**
-   * Call the actual INCL routine.
-   */
-
+  //  G4cout <<"Before PNU:" << G4endl;
   //  randomGenerator->printSeeds();
-  G4double jremx, jremy, jremz;
-  pnu(&ibert, &nopart,&izrem,&iarem,&esrem,&erecrem,&alrem,&berem,&garem,&bimpac,&jrem,
-      &jremx, &jremy, &jremz);
+  // Call the actual INCL routine:
+  pnu(&ibert, &nopart,&izrem,&iarem,&esrem,&erecrem,&alrem,&berem,&garem,&bimpac,
+      &jrem, &xjrem, &yjrem, &zjrem);
+  //  G4cout <<"After PNU:" << G4endl;
+  //  randomGenerator->printSeeds();
+  G4double mrem = int(zjrem/197.328); // CHECK
+  if (mrem > jrem) mrem=jrem;
+  if (mrem < -jrem) mrem=-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
-  G4double jprf = 0.0;           // angular momentum of the prefragment
-
-  // Mean angular momentum of prefragment.                                  
-  jprf = 0.165 * std::pow(at,(2.0/3.0)) * aprf * (at - aprf) / (at - 1.0);                               
+  G4double aprf = double(iarem);    // mass number of the prefragment
+  G4double jprf = 0.0;                // angular momentum of the prefragment
+
+  // Mean angular momentum of prefragment                                  
+  jprf = 0.165 * std::pow(at,(2.0/3.0)) * aprf*(at - aprf)/(at - 1.0);                               
   if (jprf < 0) {
     jprf = 0.0;
   }
 
-  // Reference M. de Jong, Ignatyuk, Schmidt Nuc. Phys A 613, p442, 7th line
+  // check m.de jong, ignatyuk, schmidt nuc.phys a 613, pg442, 7th line
   jprf = std::sqrt(2*jprf);
+
   jprf = jrem;
-  varntp->jremn = jrem; // Copy jrem to output tuple.
-
-  G4double numpi = 0;  // Number  of pions.
-  G4double multn = 0;  // Number (multiplicity) of neutrons.
-  G4double multp = 0;  // Number (multiplicity) of protons.
-
-  // Ecriture dans le ntuple des particules de cascade (sauf remnant).      
-  varntp->ntrack = nopart; // Nombre de particules pour ce tir.
+  varntp->jremn = jrem;      // jrem copie dans le ntuple
+
+  G4double numpi = 0;  // compteurs de pions, neutrons protons
+  G4double multn = 0; 
+  G4double multp = 0;
+
+  // ecriture dans le ntuple des particules de cascade (sauf remnant)      
+  varntp->ntrack = nopart;          // nombre de particules pour ce tir
+  if(varntp->ntrack >= VARNTPSIZE) {
+    if(verboseLevel > 2) {
+      G4cout <<"G4Incl error: Output data structure not big enough." << G4endl;
+    }
+  }
   varntp->massini = iarem;
   varntp->mzini = izrem;
@@ -533,14 +680,10 @@
   varntp->bimpact = bimpac;
   
-  /**
-   * Three ways to compute the mass of the remnant: 
-   * -from the output of the cascade and the canonic mass
-   * -from energy balance (input - all emitted energies)
-   * -following the approximations of the cugnon code (esrem...)
-   */
-  G4double f0 = calincl->targetA();
-  G4double f1 = calincl->targetZ();
-  G4double f2 = calincl->bulletE();
-  G4double mcorem = mprojo + f2 + abla->pace2(f0, f1) + f0 * uma - f1 * melec;
+  //  three ways to compute the mass of the remnant: 
+  //                -from the output of the cascade and the canonic mass
+  //                -from energy balance (input - all emitted energies)
+  //                -following the approximations of the cugnon code (esrem...)
+  G4double mcorem = mprojo + calincl->bulletE() + abla->pace2(double(calincl->targetA()),double(calincl->targetZ()))
+    + calincl->targetA()*uma - calincl->targetZ()*melec;
 
   G4double pxbil = 0.0;
@@ -549,11 +692,40 @@
 
   if(nopart > -1) {
-    for(G4int j = 0; j < nopart; j++) {
-      varntp->itypcasc[j] = 1;
+    // Fill the projectile spectator variables
+    varntp->masp = ps->a_projspec;
+    varntp->mzsp = ps->z_projspec;
+    varntp->exsp = ps->ex_projspec;
+    varntp->spectatorP1 = ps->p1_projspec;
+    varntp->spectatorP2 = ps->p2_projspec;
+    varntp->spectatorP3 = ps->p3_projspec;
+    varntp->spectatorT = ps->t_projspec;
+    for(G4int j = 0; j <= nopart; j++) {
+      if(ep[j] < 0.0) { // Workaround to avoid negative energies (and taking std::sqrt of a negative number).
+	G4cout <<"G4Incl: Not registering particle with energy: " << ep[j] << G4endl;
+	continue;
+      }
+      if(kind[j] == 0) continue; // Empty particle rows are sometimes produced by lurking indexing problems. We can simply skip those "bad" entries...
+      if(gam[j] > CLHEP::pi) {
+	if(verboseLevel > 2) {
+	  G4cout <<"G4Incl: Just avoided floating point exception by using an ugly hack..." << G4endl;
+	}
+	continue; // Avoid floating point exception
+      }
+
+      varntp->itypcasc[j] = 1; // Particle was produced by the cascade
+      // Spectators of composite projectiles (7/2006, AB)
+      // (kind is negative in that case)
+      if(kind[j] <= 0) { // Particle is a projectile spectator that comes directly from the cascade
+	kind[j] *= -1;
+	varntp->itypcasc[j]=-1;
+	//	G4cout <<"Spectator registered!" << G4endl;
+	//	continue;
+      }
+	
       // kind(): 1=proton, 2=neutron, 3=pi+, 4=pi0, 5=pi -      
       if(kind[j] == 1) { 
 	varntp->avv[j] = 1;
 	varntp->zvv[j] = 1;
-	varntp->plab[j] = std::sqrt(ep[j] * (ep[j] + 1876.5592)); // cugnon
+	varntp->plab[j] = std::sqrt(ep[j]*(ep[j]+1876.5592)); // cugnon
 	multp = multp + 1;
 	mcorem = mcorem - ep[j] - 938.27231;
@@ -568,5 +740,4 @@
 	varntp->zvv[j] = 0;
 	varntp->plab[j] = std::sqrt(ep[j]*(ep[j]+1876.5592)); // cugnon
-	//varntp->plab[j] = std::sqrt(ep[j] * (ep[j] + 1879.13126)); // PK mass check
 	multn = multn + 1;
 	mcorem = mcorem - ep[j] - 939.56563;
@@ -576,5 +747,5 @@
 	}
       }
-
+      
       if(kind[j] == 3) { 
 	varntp->avv[j] = -1;
@@ -656,5 +827,5 @@
 	mcorem = mcorem - ep[j] - 3724.818;
 	if(verboseLevel > 3) {
-	  G4cout <<"G4Incl: He4 produced! " << G4endl;
+	  G4cout <<"G4Incl: He3 produced! " << G4endl;
 	  G4cout <<"G4Incl: Momentum: "<< varntp->plab[j] << G4endl;
 	}
@@ -669,26 +840,36 @@
 
       if(verboseLevel > 3) {
-	G4cout <<"Momentum: " << varntp->plab[j] << G4endl;
-	G4cout <<"Theta: " << varntp->tetlab[j] << G4endl;
-	G4cout <<"Phi: " << varntp->philab[j] << G4endl;
+	G4cout <<"Momentum: " << varntp->plab[j]   << G4endl;
+	G4cout <<"Theta: "    << varntp->tetlab[j] << G4endl;
+	G4cout <<"Phi: "      << varntp->philab[j] << G4endl;
       }
     }
 
     // calcul de la masse (impulsion) du remnant coherente avec la conservation d'energie:
-    pcorem=std::sqrt(erecrem*(erecrem +2.*938.2796*iarem));   // cugnon
-    mcorem = 938.2796*iarem;                               // cugnon
-                
+    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
     // masse MCOREM comme ci-dessus (fausse de ~1GeV par rapport aux tables...)        
-    pxrem=pcorem*alrem;
-    pyrem=pcorem*berem;
-    pzrem=pcorem*garem;
-        
-    pxbil=pxbil+pxrem;
-    pybil=pybil+pyrem;
-    pzbil=pzbil+pzrem;
-
-    if((std::fabs(pzbil-pbeam) > 5.0) || (std::sqrt(std::pow(pxbil,2)+std::pow(pybil,2)) >= 3.0)) {
+    pxrem = pcorem * alrem;
+    pyrem = pcorem * berem;
+    pzrem = pcorem * garem;
+    varntp->pxrem = pxrem;
+    varntp->pyrem = pyrem;
+    varntp->pzrem = pzrem;
+    pxbil = pxbil + pxrem;
+    pybil = pybil + pyrem;
+    pzbil = pzbil + pzrem;
+
+    // If on purpose, add here the spectator nuclei:	
+    if(calincl->bulletType() < 0 && ps->a_projspec != 0) {
+      pxbil=pxbil+ps->p1_projspec;
+      pybil=pybil+ps->p2_projspec;
+      pzbil=pzbil+ps->p3_projspec;
+    }
+
+    if((std::fabs(pzbil - pbeam) > 5.0) || (std::sqrt(std::pow(pxbil,2) + std::pow(pybil,2)) >= 3.0)) {
       if(verboseLevel > 3) {
 	G4cout <<"Bad momentum conservation after INCL:" << G4endl;
@@ -704,35 +885,35 @@
     varntp->iafis = 0;
 
-    //  varntp->ntrack = varntp->ntrack + 1;  // on recopie le remnant dans le ntuple
+    // on recopie le remnant dans le ntuple
+    varntp->ntrack = varntp->ntrack + 1;
     varntp->massini = iarem;
     varntp->mzini = izrem;
     varntp->exini = esrem;
-    varntp->itypcasc[varntp->ntrack] = 1;
-    varntp->avv[varntp->ntrack] = iarem;
-    varntp->zvv[varntp->ntrack]= izrem;
-    varntp->plab[varntp->ntrack] = pcorem;
-    varntp->enerj[varntp->ntrack] = std::sqrt(std::pow(pcorem,2) + std::pow(mcorem,2)) - mcorem;
-    varntp->tetlab[varntp->ntrack] = 180.0*std::acos(garem)/3.141592654;
-    varntp->philab[varntp->ntrack] = 180.0*std::atan2(berem,alrem)/3.141592654;
-    varntp->ntrack++;
-    varntp->mulncasc = varntp->ntrack;
-    varntp->mulnevap = 0;
-    varntp->mulntot = varntp->mulncasc + varntp->mulnevap;
-    if(verboseLevel > 3) {
-      G4cout <<"G4Incl: Returning nucleus fragment. " << G4endl;
-      G4cout <<"G4Incl: Fragment A = " << varntp->avv[varntp->ntrack] << " Z = " << varntp->zvv[varntp->ntrack] << G4endl;
-      G4cout <<"Energy: " << varntp->enerj[varntp->ntrack] << G4endl;
-      G4cout <<"Momentum: " << varntp->plab[varntp->ntrack] << G4endl;
-      G4cout <<"Theta: " << varntp->tetlab[varntp->ntrack] << G4endl;
-      G4cout <<"Phi: " << varntp->philab[varntp->ntrack] << G4endl;
-    }
-  }
-  else {
-    if(nopart == -2) {
-      varntp->ntrack = -2; //FIX: Error flag to remove events containing unphysical events (Ekin > Ebullet).
-    }
-    else {
-      varntp->ntrack = -1;
-    }
+    varntp->pxrem = pxrem;
+    varntp->pyrem = pyrem;
+    varntp->pzrem = pzrem;
+    varntp->mcorem = mcorem;
+    varntp->erecrem = pcorem;
+    varntp->erecrem = erecrem;
+
+#ifdef G4INCLDEBUG
+    theLogger->fillHistogram1D("bimpact", varntp->bimpact);
+#endif
+
+#ifdef G4INCLDEBUG
+    theLogger->fillHistogram1D("mzini", varntp->mzini);
+#endif
+  }
+  if(nopart == -2) {
+    varntp->ntrack = -2; //FIX: Error flag to remove events containing unphysical events (Ekin > Ebullet).
+    evaporationResult->ntrack = -2; //FIX: Error flag to remove events containing unphysical events (Ekin > Ebullet).
+  }
+  else if(nopart == -1) {
+    varntp->ntrack = -1;
+    evaporationResult->ntrack = -1;
+  }
+  if(verboseLevel > 2) {
+    G4cout << __FILE__ << ":" << __LINE__ << "Dump varntp after combining: " << G4endl;
+    varntp->dump();
   }
 }
@@ -8150,7 +8331,6 @@
 }
 
-void G4Incl::standardRandom(G4double *rndm, G4long *seed)
-{
-  (*seed) = (*seed); // Avoid warning during compilation.
+void G4Incl::standardRandom(G4double *rndm, G4long*)
+{
   // Use Geant4 G4UniformRand
   //  (*rndm) = G4UniformRand();
@@ -8550,5 +8730,5 @@
     } // enddo
     G4double p_spec2=std::pow(p1_spec,2)+std::pow(p2_spec,2)+std::pow(p3_spec,2);
-    G4double s_spec = sqrt(std::pow(e_spec,2)-p_spec2);
+    G4double s_spec = std::sqrt(std::pow(e_spec,2)-p_spec2);
 
     // no projectile spectator if a>=4 and a=z or a=n (no dresner breakup)
@@ -8844,4 +9024,17 @@
 }
 
+// Initialization helper
+void G4Incl::clearState() {
+  G4int epSize = 300;
+  for(G4int i = 0; i < epSize; ++i) {
+    kind[i] = 0;
+    ep[i] = 0.0;
+    alpha[i] = 0.0;
+    beta[i] = 0.0;
+    gam[i] = 0.0;
+  }
+  inside_step = 0;
+}
+
 // C------------------------------------------------------------------------
 // C Logging
@@ -8923,4 +9116,9 @@
 //       COMMON/BL2/CROIS(19900),K,IND(20000),JND(20000) 
 //       COMMON/BL3/R1,R2,X1(300),X2(300),X3(300),IA1,IA2,RAB2             P-N22270
+
+  if(index < 0) {
+    G4cout <<";; Error! index < 0! index = " << index << G4endl;
+    return;
+  }
 
   G4int i_ind1 = bl1->ind1[bl2->ind[index]];
Index: trunk/source/processes/hadronic/models/incl/src/G4InclAblaCascadeInterface.cc
===================================================================
--- trunk/source/processes/hadronic/models/incl/src/G4InclAblaCascadeInterface.cc	(revision 1340)
+++ trunk/source/processes/hadronic/models/incl/src/G4InclAblaCascadeInterface.cc	(revision 1347)
@@ -24,5 +24,5 @@
 // ********************************************************************
 //
-// $Id: G4InclAblaCascadeInterface.cc,v 1.16 2010/10/29 06:48:43 gunter Exp $ 
+// $Id: G4InclAblaCascadeInterface.cc,v 1.20 2010/11/17 20:19:09 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -116,6 +116,6 @@
   incl->setInput(calincl);
 
-  G4InclInput::printProjectileTargetInfo(aTrack, theNucleus);
-  calincl->printInfo();
+  //  G4InclInput::printProjectileTargetInfo(aTrack, theNucleus);
+  //  calincl->printInfo();
 
 #ifdef DEBUGINCL
@@ -439,4 +439,7 @@
 	}
       }
+      delete theFermiBreakupResult;
+      theFermiBreakupResult = 0;
+
       if(std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 
 	G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl;
Index: trunk/source/processes/hadronic/models/incl/src/G4InclAblaDataFile.cc
===================================================================
--- trunk/source/processes/hadronic/models/incl/src/G4InclAblaDataFile.cc	(revision 1340)
+++ trunk/source/processes/hadronic/models/incl/src/G4InclAblaDataFile.cc	(revision 1347)
@@ -24,5 +24,5 @@
 // ********************************************************************
 //
-// $Id: G4InclAblaDataFile.cc,v 1.9 2010/06/14 16:10:01 gcosmo Exp $ 
+// $Id: G4InclAblaDataFile.cc,v 1.10 2010/11/17 20:19:09 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -116,5 +116,4 @@
   vgsldin.close();
 
-  int A = 0, Zbegin = 0, Zend = 0;
   G4String str1, str2, str3;
   for(int i = 0; i < 500; i++) {
@@ -124,9 +123,13 @@
   }
   
+  int A = 0, Zbegin = 0, Zend = 0;
   for(int i = 0; i < massnumbers; i++) {
     pace2in >> str1 >> A >> str2 >> Zbegin >> str3 >> Zend;
-    for(int j = Zbegin; j <= Zend; j++) {
-      pace2in >> pace2;
-      setPace2(A, j, pace2);
+    if(Zbegin >= 0 && Zbegin < getPaceCols() &&
+       A >= 0 && A < getPaceRows()) {
+      for(int j = Zbegin; j <= Zend; j++) {
+	pace2in >> pace2;
+	setPace2(A, j, pace2);
+      }
     }
   }
Index: trunk/source/processes/hadronic/models/incl/src/G4InclAblaLightIonInterface.cc
===================================================================
--- trunk/source/processes/hadronic/models/incl/src/G4InclAblaLightIonInterface.cc	(revision 1340)
+++ trunk/source/processes/hadronic/models/incl/src/G4InclAblaLightIonInterface.cc	(revision 1347)
@@ -24,5 +24,5 @@
 // ********************************************************************
 //
-// $Id: G4InclAblaLightIonInterface.cc,v 1.14 2010/10/29 06:48:43 gunter Exp $ 
+// $Id: G4InclAblaLightIonInterface.cc,v 1.16 2010/11/17 20:19:09 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -103,21 +103,18 @@
   }
 
-  if(verboseLevel > 1) {
-    G4Calincl::printProjectileTargetInfo(aTrack, theNucleus);
-  }
-
   // Inverse kinematics for targets with Z = 1 and A = 1
   //  if(false) {
   G4LorentzRotation toBreit = aTrack.Get4Momentum().boostVector();
 
-  if(theNucleus.GetZ_asInt() == 1 && theNucleus.GetA_asInt() == 1 && G4Calincl::canUseInverseKinematics(aTrack, theNucleus)) {
+  if(theNucleus.GetZ_asInt() == 1 && theNucleus.GetA_asInt() == 1 && G4InclInput::canUseInverseKinematics(aTrack, theNucleus)) {
     G4ParticleDefinition *oldTargetDef = theTableOfParticles->GetIon(theNucleus.GetA_asInt(), theNucleus.GetZ_asInt(), 0.0);
     const G4ParticleDefinition *oldProjectileDef = aTrack.GetDefinition();
 
+    if(oldProjectileDef != 0 && oldTargetDef != 0) {
     G4int oldTargetA = oldTargetDef->GetAtomicMass();
     G4int newTargetA = oldProjectileDef->GetAtomicMass();
     G4int newTargetZ = oldProjectileDef->GetAtomicNumber();
 
-    if(newTargetA > 0 && newTargetZ > 0 && oldTargetDef != 0 && oldProjectileDef != 0) {
+    if(newTargetA > 0 && newTargetZ > 0) {
       G4Nucleus swappedTarget(oldProjectileDef->GetAtomicMass(), oldProjectileDef->GetAtomicNumber());
 
@@ -137,4 +134,5 @@
       G4cout <<"Badly defined target after swapping. Falling back to normal (non-swapped) mode." << G4endl;
       calincl = new G4InclInput(aTrack, theNucleus, false);
+    }
     }
   } else {
@@ -489,4 +487,7 @@
 	}
       }
+      delete theSpectatorFermiBreakupResult;
+      theSpectatorFermiBreakupResult = 0;
+
       if(std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 
 	G4cout <<"Four-momentum balance after spectator nucleus Fermi break-up:" << G4endl;
@@ -588,4 +589,7 @@
 	}
       }
+      delete theFermiBreakupResult;
+      theFermiBreakupResult = 0;
+
       if(std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 
 	G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl;
@@ -697,4 +701,5 @@
   }
 
+  delete fermiBreakUp;
   delete calincl;
   calincl = 0;
Index: trunk/source/processes/hadronic/models/incl/src/G4InclCascadeInterface.cc
===================================================================
--- trunk/source/processes/hadronic/models/incl/src/G4InclCascadeInterface.cc	(revision 1340)
+++ trunk/source/processes/hadronic/models/incl/src/G4InclCascadeInterface.cc	(revision 1347)
@@ -24,5 +24,5 @@
 // ********************************************************************
 //
-// $Id: G4InclCascadeInterface.cc,v 1.12 2010/10/26 02:47:59 kaitanie Exp $ 
+// $Id: G4InclCascadeInterface.cc,v 1.15 2010/11/17 20:19:09 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -34,9 +34,11 @@
 
 #include "G4InclCascadeInterface.hh"
+#include "G4FermiBreakUp.hh"
 #include "math.h"
 #include "G4GenericIon.hh"
 #include "CLHEP/Random/Random.h"
 
-G4InclCascadeInterface::G4InclCascadeInterface()
+G4InclCascadeInterface::G4InclCascadeInterface(const G4String& nam)
+  :G4VIntraNuclearTransportModel(nam)
 {
   hazard = new G4Hazard();
@@ -45,8 +47,15 @@
 
   varntp = new G4VarNtp();
-  inclInput = 0;
+  calincl = 0;
   ws = new G4Ws();
   mat = new G4Mat();
-  incl = new G4Incl(hazard, inclInput, ws, mat, varntp);
+  incl = new G4Incl(hazard, calincl, ws, mat, varntp);
+
+  theExcitationHandler = new G4ExcitationHandler;
+  thePrecoModel = new G4PreCompoundModel(theExcitationHandler);
+
+  if(!getenv("G4INCLABLANOFERMIBREAKUP")) { // Use Fermi Break-up by default if it is NOT explicitly disabled
+    incl->setUseFermiBreakUp(true);
+  }
 
   verboseLevel = 0;
@@ -55,7 +64,9 @@
 G4InclCascadeInterface::~G4InclCascadeInterface()
 {
+  delete thePrecoModel;
+  delete theExcitationHandler;
+
   delete hazard;
   delete varntp;
-  delete inclInput;
   delete ws;
   delete mat;
@@ -68,4 +79,5 @@
 
   G4int particleI;
+  G4int bulletType = 0;
 
   // Print diagnostic messages: 0 = silent, 1 and 2 = verbose
@@ -91,4 +103,7 @@
   G4DynamicParticle *cascadeParticle = 0;
   G4ParticleDefinition *aParticleDefinition = 0;
+  
+  G4ReactionProductVector *thePrecoResult = 0;
+  G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
 
   // INCL assumes the projectile particle is going in the direction of
@@ -102,7 +117,11 @@
   G4LorentzRotation toLabFrame = toZ.inverse();
 
-  inclInput = new G4InclInput(aTrack, theNucleus, false);
-
   theResult.Clear(); // Make sure the output data structure is clean.
+
+  calincl = new G4InclInput(aTrack, theNucleus, false);
+  incl->setInput(calincl);
+
+  //  G4InclInput::printProjectileTargetInfo(aTrack, theNucleus);
+  //  calincl->printInfo();
 
 #ifdef DEBUGINCL
@@ -117,9 +136,10 @@
   G4double eKinSum = bulletE;
   G4LorentzVector labv = G4LorentzVector(0.0, 0.0, std::sqrt(bulletE*(bulletE + 2.*mass)), bulletE + mass + amass);
+  G4LorentzVector labvA = G4LorentzVector(0.0, 0.0, 0.0, 0.0);
   G4cout <<"Energy in the beginning = " << labv.e() / MeV << G4endl;
 #endif
 
   // Check wheter the input is acceptable.
-  if((inclInput->bulletType() != 0) && ((inclInput->targetA() != 1) && (inclInput->targetZ() != 1))) {
+  if((calincl->bulletType() != 0) && ((calincl->targetA() != 1) && (calincl->targetZ() != 1))) {
     ws->nosurf = -2;  // Nucleus surface, -2 = Woods-Saxon 
     ws->xfoisa = 8;
@@ -130,6 +150,6 @@
 
     mat->nbmat = 1;
-    mat->amat[0] = int(inclInput->targetA());
-    mat->zmat[0] = int(inclInput->targetZ());
+    mat->amat[0] = int(calincl->targetA());
+    mat->zmat[0] = int(calincl->targetZ());
 
     incl->initIncl(true);
@@ -140,5 +160,5 @@
         G4cout <<"G4InclCascadeInterface: Try number = " << nTries << G4endl; 
       }
-      incl->processEventIncl(inclInput);
+      incl->processEventIncl(calincl);
 
       if(verboseLevel > 1) {
@@ -151,38 +171,18 @@
        * Diagnostic output
        */
-      G4cout <<"G4InclCascadeInterface: Bullet type: " << inclInput->bulletType() << G4endl;
-      G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << inclInput->bulletE() << " MeV" << G4endl;
-
-      G4cout <<"G4InclCascadeInterface: Target A:  " << inclInput->targetA() << G4endl;
-      G4cout <<"G4InclCascadeInterface: Target Z:  " << inclInput->targetZ() << G4endl;
+      G4cout <<"G4InclCascadeInterface: Bullet type: " << calincl->bulletType() << G4endl;
+      G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl;
+
+      G4cout <<"G4InclCascadeInterface: Target A:  " << calincl->targetA() << G4endl;
+      G4cout <<"G4InclCascadeInterface: Target Z:  " << calincl->targetZ() << G4endl;
 
       if(verboseLevel > 3) {
-        diagdata <<"G4InclCascadeInterface: Bullet type: " << inclInput->bulletType() << G4endl;
-        diagdata <<"G4InclCascadeInterface: Bullet energy: " << inclInput->bulletE() << " MeV" << G4endl;
+        diagdata <<"G4InclCascadeInterface: Bullet type: " << calincl->bulletType() << G4endl;
+        diagdata <<"G4InclCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl;
         
-        diagdata <<"G4InclCascadeInterface: Target A:  " << inclInput->targetA() << G4endl;
-        diagdata <<"G4InclCascadeInterface: Target Z:  " << inclInput->targetZ() << G4endl;
-      }
-
-      // for(particleI = 0; particleI < varntp->ntrack; particleI++) {
-      //   G4cout << n                       << " " << inclInput->f[6]             << " " << inclInput->f[2] << " ";
-      //   G4cout << varntp->massini         << " " << varntp->mzini             << " ";
-      //   G4cout << varntp->exini           << " " << varntp->mulncasc          << " " << varntp->mulnevap          << " " << varntp->mulntot << " ";
-      //   G4cout << varntp->bimpact         << " " << varntp->jremn             << " " << varntp->kfis              << " " << varntp->estfis << " ";
-      //   G4cout << varntp->izfis           << " " << varntp->iafis             << " " << varntp->ntrack            << " " << varntp->itypcasc[particleI] << " ";
-      //   G4cout << varntp->avv[particleI]  << " " << varntp->zvv[particleI]    << " " << varntp->enerj[particleI]  << " ";
-      //   G4cout << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl;
-      //   // For diagnostic output
-      //   if(verboseLevel > 3) {
-      //     diagdata << n                       << " " << inclInput->f[6]             << " " << inclInput->f[2] << " ";
-      //     diagdata << varntp->massini         << " " << varntp->mzini             << " ";
-      //     diagdata << varntp->exini           << " " << varntp->mulncasc          << " " << varntp->mulnevap          << " " << varntp->mulntot << " ";
-      //     diagdata << varntp->bimpact         << " " << varntp->jremn             << " " << varntp->kfis              << " " << varntp->estfis << " ";
-      //     diagdata << varntp->izfis           << " " << varntp->iafis             << " " << varntp->ntrack            << " ";
-      // 	  diagdata                                                                       << varntp->itypcasc[particleI] << " ";
-      //     diagdata << varntp->avv[particleI]  << " " << varntp->zvv[particleI]    << " " << varntp->enerj[particleI]  << " ";
-      //     diagdata << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl;
-      //   }
-      // }
+        diagdata <<"G4InclCascadeInterface: Target A:  " << calincl->targetA() << G4endl;
+        diagdata <<"G4InclCascadeInterface: Target Z:  " << calincl->targetZ() << G4endl;
+      }
+
     }
 
@@ -196,10 +196,14 @@
 
       theResult.SetStatusChange(stopAndKill);
-      aParticleDefinition = G4InclInput::getParticleDefinition(inclInput->bulletType());
-
-      cascadeParticle = new G4DynamicParticle();
-      cascadeParticle->SetDefinition(aParticleDefinition);
-      cascadeParticle->Set4Momentum(aTrack.Get4Momentum());
-      theResult.AddSecondary(cascadeParticle);
+      
+      G4int bulletType = calincl->bulletType();
+      aParticleDefinition = G4InclInput::getParticleDefinition(bulletType);
+
+      if(aParticleDefinition != 0) {
+	cascadeParticle = new G4DynamicParticle();
+	cascadeParticle->SetDefinition(aParticleDefinition);
+	cascadeParticle->Set4Momentum(aTrack.Get4Momentum());
+	theResult.AddSecondary(cascadeParticle);
+      }
     }
 
@@ -209,5 +213,12 @@
 
 #ifdef DEBUGINCL
-    G4cout << "E [MeV]" << std::setw(12) << " Ekin [MeV]" << std::setw(12) << " E* [MeV]" << std::setw(12) << "Px [MeV]" << std::setw(12) << " Py [MeV]" << std::setw(12) << "Pz [MeV]" << std::setw(12) << "Pt [MeV]" << std::setw(12) << "A" << std::setw(12) << "Z" << G4endl;  
+    G4cout << "E [MeV]" << std::setw(12)
+	   << " Ekin [MeV]" << std::setw(12)
+	   << "Px [MeV]" << std::setw(12)
+	   << " Py [MeV]" << std::setw(12)
+	   << "Pz [MeV]" << std::setw(12)
+	   << "Pt [MeV]" << std::setw(12)
+	   << "A" << std::setw(12)
+	   << "Z" << G4endl;
 #endif
 
@@ -276,5 +287,4 @@
       if((varntp->avv[particleI] > 1) && (varntp->zvv[particleI] >= 1)) { // Nucleus fragment
         G4ParticleDefinition * aIonDef = 0;
-        G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
 
         G4int A = G4int(varntp->avv[particleI]);
@@ -315,7 +325,10 @@
 	G4double p = mom.mag();
 	labv -= fm;
- 	G4double px = mom.x() * MeV;
- 	G4double py = mom.y() * MeV;
- 	G4double pz = mom.z() * MeV;
+	if(varntp->avv[particleI] > 1) {
+	  labvA += fm;
+	}
+	G4double px = mom.x() * MeV;
+	G4double py = mom.y() * MeV;
+	G4double pz = mom.z() * MeV;
 	G4double pt = std::sqrt(px*px+py*py);
 	G4double e  = fm.e();
@@ -330,5 +343,4 @@
 	G4cout << fm.e() / MeV
 	       << std::setw(12) << cascadeParticle->GetKineticEnergy() / MeV
-	       << std::setw(12) << exE / MeV
 	       << std::setw(12) << mom.x() / MeV
 	       << std::setw(12) << mom.y() / MeV
@@ -351,10 +363,113 @@
       }
     }
+
+      G4double nuclearMass = G4NucleiProperties::GetNuclearMass(G4int(varntp->massini), G4int(varntp->mzini)) + varntp->exini * MeV;
+      G4LorentzVector fragmentMomentum(varntp->pxrem * MeV, varntp->pyrem * MeV, varntp->pzrem * MeV,
+				       varntp->erecrem * MeV + nuclearMass);
+      G4double momentumScaling = G4InclUtils::calculate4MomentumScaling(G4int(varntp->massini), G4int(varntp->mzini),
+									varntp->exini,
+									varntp->erecrem,
+									varntp->pxrem,
+									varntp->pyrem,
+									varntp->pzrem);
+      G4LorentzVector p4(momentumScaling * varntp->pxrem * MeV, momentumScaling * varntp->pyrem * MeV,
+			 momentumScaling * varntp->pzrem * MeV,
+      			 varntp->erecrem + nuclearMass);
+
+      // For four-momentum, baryon number and charge conservation check:
+      G4LorentzVector fourMomentumBalance = p4;
+      G4int baryonNumberBalance = G4int(varntp->massini);
+      G4int chargeBalance = G4int(varntp->mzini);
+
+      G4LorentzRotation toFragmentZ;
+      toFragmentZ.rotateZ(-p4.theta());
+      toFragmentZ.rotateY(-p4.phi());
+      G4LorentzRotation toFragmentLab = toFragmentZ.inverse();
+      p4 *= toFragmentZ;
+
+      G4LorentzVector p4rest = p4;
+      p4rest.boost(-p4.boostVector());
+      if(verboseLevel > 0) {
+	G4cout <<"Cascade remnant nucleus:" << G4endl;
+	G4cout <<"p4: " << G4endl;
+	G4cout <<" px: " << p4.px() <<" py: " << p4.py() <<" pz: " << p4.pz() << G4endl;
+	G4cout <<" E = " << p4.e() << G4endl;
+
+	G4cout <<"p4rest: " << G4endl;
+	G4cout <<" px: " << p4rest.px() <<" py: " << p4rest.py() <<" pz: " << p4rest.pz() << G4endl;
+	G4cout <<" E = " << p4rest.e() << G4endl;
+      }
+
+      G4Fragment theCascadeRemnant(G4int(varntp->massini), G4int(varntp->mzini), p4rest);
+      thePrecoResult = thePrecoModel->DeExcite(theCascadeRemnant);
+      if(thePrecoResult != 0) {
+      G4ReactionProductVector::iterator fragment;
+      for(fragment = thePrecoResult->begin(); fragment != thePrecoResult->end(); fragment++) {
+	G4ParticleDefinition *theFragmentDefinition = (*fragment)->GetDefinition();
+
+	if(theFragmentDefinition != 0) {
+	  G4DynamicParticle *theFragment = new G4DynamicParticle(theFragmentDefinition, (*fragment)->GetMomentum());
+	  G4LorentzVector labMomentum = theFragment->Get4Momentum();
+	  labMomentum.boost(p4.boostVector());
+	  labMomentum *= toFragmentLab;
+	  labMomentum *= toLabFrame;
+	  theFragment->Set4Momentum(labMomentum);
+	  fourMomentumBalance -= theFragment->Get4Momentum();
+	  baryonNumberBalance -= theFragmentDefinition->GetAtomicMass();
+	  chargeBalance -= theFragmentDefinition->GetAtomicNumber();
+	  if(verboseLevel > 0) {
+	    G4cout <<"Resulting fragment: " << G4endl;
+	    G4cout <<" kinetic energy = " << theFragment->GetKineticEnergy() / MeV << " MeV" << G4endl;
+	    G4cout <<" momentum = " << theFragment->GetMomentum().mag() / MeV << " MeV" << G4endl;
+	  }
+	  theResult.AddSecondary(theFragment);
+	} else {
+	  G4cout <<"G4InclCascadeInterface: Error. Fragment produced by Fermi break-up does not exist." << G4endl;
+	  G4cout <<"Resulting fragment: " << G4endl;
+	  G4cout <<" momentum = " << (*fragment)->GetMomentum().mag() / MeV << " MeV" << G4endl;
+	}
+      }
+      delete thePrecoResult;
+      thePrecoResult = 0;
+
+      if(verboseLevel > 1 && std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 
+	G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl;
+	G4cout <<"Magnitude: " << fourMomentumBalance.mag() / MeV << " MeV" << G4endl;
+	G4cout <<"Vector components (px, py, pz, E) = ("
+	       << fourMomentumBalance.px() << ", "
+	       << fourMomentumBalance.py() << ", "
+	       << fourMomentumBalance.pz() << ", "
+	       << fourMomentumBalance.e() << ")" << G4endl;
+      }
+      if(baryonNumberBalance != 0 && verboseLevel > 1) {
+	G4cout <<"Baryon number balance after remnant nucleus Fermi break-up: " << baryonNumberBalance << G4endl;
+      }
+      if(chargeBalance != 0 && verboseLevel > 1) {
+	G4cout <<"Charge balance after remnant nucleus Fermi break-up: " << chargeBalance << G4endl;
+      }
+      }
+      //    } // if(needsFermiBreakUp)
+
 #ifdef DEBUGINCL
     G4cout <<"--------------------------------------------------------------------------------" << G4endl;
     G4double pt = std::sqrt(std::pow(labv.x(), 2) + std::pow(labv.y(), 2));
-    G4cout << labv.e() / MeV << std::setw(12) << eKinSum / MeV << std::setw(12) << labv.x() << std::setw(12) << labv.y() << std::setw(12) << labv.z() << std::setw(12) <<  pt / MeV << std::setw(12) << baryonNumber << std::setw(12) << chargeNumber << " totals" << G4endl;
+    G4double ptA = std::sqrt(std::pow(labvA.x(), 2) + std::pow(labvA.y(), 2));
+    G4cout << labv.e() / MeV << std::setw(12)
+	   << eKinSum  / MeV << std::setw(12)
+	   << labv.x() / MeV << std::setw(12)
+	   << labv.y() / MeV << std::setw(12)
+	   << labv.z() / MeV << std::setw(12)
+	   << pt       / MeV << std::setw(12)
+	   << baryonNumber << std::setw(12)
+	   << chargeNumber << " totals" << G4endl;
+    G4cout << " - " << std::setw(12)
+	   << " - " << std::setw(12)
+	   << labvA.x() / MeV << std::setw(12)
+	   << labvA.y() / MeV << std::setw(12)
+	   << labvA.z() / MeV << std::setw(12)
+	   << ptA       / MeV << std::setw(12)
+	   << " - " << std::setw(12) << " - " << " totals ABLA" << G4endl;
     G4cout << G4endl;
-
+ 
     if(verboseLevel > 3) {
       if(baryonNumber != 0) {
@@ -369,5 +484,5 @@
     }
 #endif
-    
+   
     varntp->ntrack = 0; // Clean up the number of generated particles in the event.
   }
@@ -379,5 +494,4 @@
     theResult.SetStatusChange(stopAndKill);
 
-    G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
     cascadeParticle = new G4DynamicParticle(theTableOfParticles->FindParticle(aTrack.GetDefinition()), aTrack.Get4Momentum());
 
@@ -388,8 +502,8 @@
     }
     if(verboseLevel > 3) {
-      diagdata <<"ERROR G4InclCascadeInterface: Processing event number (internal) failed " << eventNumber << G4endl;
-    }
-
-    if(inclInput->bulletType() == 0) {
+      diagdata <<"ERROR G4InclCascadeInterface: Error processing event number (internal) failed " << eventNumber << G4endl;
+    }
+
+    if(bulletType == 0) {
       if(verboseLevel > 1) {
 	G4cout <<"G4InclCascadeInterface: Unknown bullet type" << G4endl;      
@@ -402,24 +516,24 @@
     }
 
-    if((inclInput->targetA() == 1) && (inclInput->targetZ() == 1)) { // Unsupported target
+    if((calincl->targetA() == 1) && (calincl->targetZ() == 1)) { // Unsupported target
       if(verboseLevel > 1) {
 	G4cout <<"Unsupported target: " << G4endl;
-	G4cout <<"Target A: " << inclInput->targetA() << G4endl;
-	G4cout <<"TargetZ: " << inclInput->targetZ() << G4endl;
+	G4cout <<"Target A: " << calincl->targetA() << G4endl;
+	G4cout <<"TargetZ: " << calincl->targetZ() << G4endl;
       }
       if(verboseLevel > 3) {
         diagdata <<"Unsupported target: " << G4endl;
-        diagdata <<"Target A: " << inclInput->targetA() << G4endl;
-	diagdata <<"TargetZ: " << inclInput->targetZ() << G4endl;
-      }
-    }
-
-    if(inclInput->bulletE() < 100) { // INCL does not support E < 100 MeV.
+        diagdata <<"Target A: " << calincl->targetA() << G4endl;
+	diagdata <<"TargetZ: " << calincl->targetZ() << G4endl;
+      }
+    }
+
+    if(calincl->bulletE() < 100) { // INCL does not support E < 100 MeV.
       if(verboseLevel > 1) {
-	G4cout <<"Unsupported bullet energy: " << inclInput->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;
+	G4cout <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;
 	G4cout <<"WARNING: Returning the original bullet with original energy back to Geant4." << G4endl;
       }
       if(verboseLevel > 3) {
-        diagdata <<"Unsupported bullet energy: " << inclInput->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;
+        diagdata <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;
       }
     }
@@ -430,4 +544,6 @@
   }
 
+  delete calincl;
+  calincl = 0;
   return &theResult;
 } 
Index: trunk/source/processes/hadronic/models/incl/src/G4InclLightIonInterface.cc
===================================================================
--- trunk/source/processes/hadronic/models/incl/src/G4InclLightIonInterface.cc	(revision 1340)
+++ trunk/source/processes/hadronic/models/incl/src/G4InclLightIonInterface.cc	(revision 1347)
@@ -24,5 +24,5 @@
 // ********************************************************************
 //
-// $Id: G4InclLightIonInterface.cc,v 1.12 2010/10/26 02:47:59 kaitanie Exp $ 
+// $Id: G4InclLightIonInterface.cc,v 1.15 2010/11/17 20:19:09 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -31,5 +31,8 @@
 // Aatos Heikkinen, HIP (project coordination)
 
+#include <vector>
+
 #include "G4InclLightIonInterface.hh"
+#include "G4FermiBreakUp.hh"
 #include "math.h"
 #include "G4GenericIon.hh"
@@ -39,6 +42,10 @@
 {
   hazard = new G4Hazard();
+
   const G4long* table_entry = CLHEP::HepRandom::getTheSeeds(); // Get random seed from CLHEP.
   hazard->ial = (*table_entry);
+
+  theExcitationHandler = new G4ExcitationHandler;
+  thePrecoModel = new G4PreCompoundModel(theExcitationHandler);
 
   varntp = new G4VarNtp();
@@ -47,10 +54,22 @@
   mat = new G4Mat();
   incl = new G4Incl(hazard, calincl, ws, mat, varntp);
-
+  useProjectileSpectator = true;
+  useFermiBreakup = true;
+  incl->setUseProjectileSpectators(useProjectileSpectator);
+  if(!getenv("G4INCLABLANOFERMIBREAKUP")) { // Use Fermi Break-up by default if it is NOT explicitly disabled
+    incl->setUseFermiBreakUp(true);
+    useFermiBreakup = true;
+  }
   verboseLevel = 0;
+  if(getenv("G4INCLVERBOSE")) {
+    verboseLevel = 1;
+  }
 }
 
 G4InclLightIonInterface::~G4InclLightIonInterface()
 {
+  delete thePrecoModel;
+  delete theExcitationHandler;
+
   delete hazard;
   delete varntp;
@@ -63,15 +82,23 @@
 G4HadFinalState* G4InclLightIonInterface::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
 {
+  //  const G4bool useFermiBreakup = false;
   G4int maxTries = 200;
 
-  G4int bulletType = 0;
   G4int particleI;
 
-  // Print diagnostic messages: 0 = silent, 1 and 2 = verbose
-  verboseLevel = 0;
+  G4int baryonNumberBalanceInINCL = 0;
+  G4int chargeNumberBalanceInINCL = 0;
+
+  G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
 
   // Increase the event number:
   eventNumber++;
 
+  // Clean up the INCL input
+  if(calincl != 0) {
+    delete calincl;
+    calincl = 0;
+  }
+
   if (verboseLevel > 1) {
     G4cout << " >>> G4InclLightIonInterface::ApplyYourself called" << G4endl;
@@ -82,9 +109,40 @@
   }
 
-  // INCL4 needs the energy in units MeV
-  G4double bulletE = aTrack.GetKineticEnergy() / MeV;
-
-  G4int targetA = theNucleus.GetA_asInt();
-  G4int targetZ = theNucleus.GetZ_asInt();
+  // Inverse kinematics for targets with Z = 1 and A = 1
+  //  if(false) {
+  G4LorentzRotation toBreit = aTrack.Get4Momentum().boostVector();
+
+  if(theNucleus.GetZ_asInt() == 1 && theNucleus.GetA_asInt() == 1 && G4InclInput::canUseInverseKinematics(aTrack, theNucleus)) {
+    G4ParticleDefinition *oldTargetDef = theTableOfParticles->GetIon(theNucleus.GetA_asInt(), theNucleus.GetZ_asInt(), 0.0);
+    const G4ParticleDefinition *oldProjectileDef = aTrack.GetDefinition();
+
+    if(oldTargetDef != 0 && oldProjectileDef != 0) {
+    G4int oldTargetA = oldTargetDef->GetAtomicMass();
+    G4int newTargetA = oldProjectileDef->GetAtomicMass();
+    G4int newTargetZ = oldProjectileDef->GetAtomicNumber();
+
+    if(newTargetA > 0 && newTargetZ > 0) {
+      G4Nucleus swappedTarget(oldProjectileDef->GetAtomicMass(), oldProjectileDef->GetAtomicNumber());
+
+      //      G4cout <<"Original projectile kinE = " << aTrack.GetKineticEnergy() / MeV << G4endl;
+
+      // We need the same energy/nucleon.
+      G4double projectileE = ((aTrack.GetKineticEnergy() / MeV) / newTargetA) * oldTargetA * MeV;
+
+      //    G4cout <<"projectileE = " << projectileE << G4endl;
+      G4DynamicParticle swappedProjectileParticle(oldTargetDef, G4ThreeVector(0.0, 0.0, 1.0), projectileE);
+      const G4LorentzVector swapped4Momentum = (swappedProjectileParticle.Get4Momentum()*=toBreit);
+      swappedProjectileParticle.Set4Momentum(swapped4Momentum);
+      const G4HadProjectile swappedProjectile(swappedProjectileParticle);
+      //  G4cout <<"New projectile kinE = " << swappedProjectile.GetKineticEnergy() / MeV << G4endl;
+      calincl = new G4InclInput(swappedProjectile, swappedTarget, true);
+    } else {
+      G4cout <<"Badly defined target after swapping. Falling back to normal (non-swapped) mode." << G4endl;
+      calincl = new G4InclInput(aTrack, theNucleus, false);
+    }
+    }
+  } else {
+    calincl = new G4InclInput(aTrack, theNucleus, false);
+  }
 
   G4double eKin;
@@ -103,20 +161,62 @@
   G4LorentzRotation toLabFrame = toZ.inverse();
 
+  /*
+  G4cout <<"Projectile theta = " << projectileMomentum.theta() << " phi = " << projectileMomentum.phi() << G4endl;
+  G4cout <<"Projectile momentum "
+	 << "(px = " << projectileMomentum.px()
+	 << ", py = " << projectileMomentum.py()
+	 << ", pz = " << projectileMomentum.pz() << ")" << G4endl;
+  G4cout << "Projectile energy = " << bulletE << " MeV" << G4endl;
+  */
+
+  G4ReactionProductVector *thePrecoResult = 0;
+  G4ReactionProductVector *theSpectatorPrecoResult = 0;
+
   theResult.Clear(); // Make sure the output data structure is clean.
+
+  std::vector<G4DynamicParticle*> result; // Temporary list for the results
 
   // Map Geant4 particle types to corresponding INCL4 types.
   enum bulletParticleType {nucleus = 0, proton = 1, neutron = 2, pionPlus = 3, pionZero = 4, 
-                           pionMinus = 5, deuteron = 6, triton = 7, he3 = 8, he4 = 9};
-
-  // Coding particles for use with INCL4 and ABLA 
-  if (aTrack.GetDefinition() == G4Deuteron::Deuteron()   ) bulletType = deuteron;
-  if (aTrack.GetDefinition() == G4Triton::Triton()       ) bulletType = triton;
-  if (aTrack.GetDefinition() == G4He3::He3()             ) bulletType = he3;
-  if (aTrack.GetDefinition() == G4Alpha::Alpha()         ) bulletType = he4;
-
-  calincl = new G4InclInput();
+                           pionMinus = 5, deuteron = 6, triton = 7, he3 = 8, he4 = 9,
+                           c12 = -12}; // Carbon beam support.
+
+  G4int bulletType = calincl->bulletType();
+  chargeNumberBalanceInINCL = calincl->targetZ();
+  baryonNumberBalanceInINCL = calincl->targetA();
+
+  //  G4cout <<"Type of the projectile (INCL projectile code): " << bulletType << G4endl;
+
+  if(bulletType == proton) {
+    chargeNumberBalanceInINCL += 1;
+    baryonNumberBalanceInINCL += 1;
+  } else if(bulletType == neutron) {
+    baryonNumberBalanceInINCL += 1;
+  } else if(bulletType == pionPlus) { //Note: positive pion doesn't contribute to the baryon and charge number counters
+    chargeNumberBalanceInINCL += 1;
+  } else if(bulletType == pionMinus) {
+    chargeNumberBalanceInINCL -= 1;
+  } else if(bulletType == deuteron) {
+    chargeNumberBalanceInINCL += 1;
+    baryonNumberBalanceInINCL += 2;
+  } else if(bulletType == triton) {
+    chargeNumberBalanceInINCL += 1;
+    baryonNumberBalanceInINCL += 3;
+  } else if(bulletType == he3) {
+    chargeNumberBalanceInINCL += 2;
+    baryonNumberBalanceInINCL += 3;
+  } else if(bulletType == he4) {
+    chargeNumberBalanceInINCL += 2;
+    baryonNumberBalanceInINCL += 4;
+  } if(bulletType == c12) {
+    chargeNumberBalanceInINCL += 6;
+    baryonNumberBalanceInINCL += 12;
+  } if(bulletType == -666) {
+    chargeNumberBalanceInINCL += calincl->extendedProjectileZ();
+    baryonNumberBalanceInINCL += calincl->extendedProjectileA();
+  }
 
   // Check wheter the input is acceptable.
-  if((bulletType != 0) && ((targetA != 1) && (targetZ != 1))) {
+  if((bulletType != 0) && ((calincl->targetA() != 1) && (calincl->targetZ() != 1))) {
     ws->nosurf = -2;  // Nucleus surface, -2 = Woods-Saxon 
     ws->xfoisa = 8;
@@ -128,6 +228,7 @@
     mat->nbmat = 1;
     mat->amat[0] = int(calincl->targetA());
-    mat->zmat[0] = int(calincl->targetZ());
-
+    mat->zmat[0] = int(calincl->targetA());
+
+    incl->setInput(calincl);
     incl->initIncl(true);
 
@@ -148,16 +249,20 @@
        * Diagnostic output
        */
-      G4cout <<"G4InclLightIonInterface: Bullet type: " << bulletType << G4endl;
-      G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << bulletE << " MeV" << G4endl;
-
-      G4cout <<"G4InclLightIonInterface: Target A:  " << targetA << G4endl;
-      G4cout <<"G4InclLightIonInterface: Target Z:  " << targetZ << G4endl;
+      G4cout <<"G4InclLightIonInterface: Bullet type: " << calincl->bulletType() << G4endl;
+      G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl;
+      if(bulletType == -666) {
+	G4cout <<"   Extended projectile: A = " << calincl->extendedProjectileA()
+	       <<" Z = " << calincl->extendedProjectileZ() << G4endl;
+      }
+
+      G4cout <<"G4InclLightIonInterface: Target A:  " << calincl->targetA() << G4endl;
+      G4cout <<"G4InclLightIonInterface: Target Z:  " << calincl->targetZ() << G4endl;
 
       if(verboseLevel > 3) {
-        diagdata <<"G4InclLightIonInterface: Bullet type: " << bulletType << G4endl;
-        diagdata <<"G4InclLightIonInterface: Bullet energy: " << bulletE << " MeV" << G4endl;
+        diagdata <<"G4InclLightIonInterface: Bullet type: " << calincl->bulletType() << G4endl;
+        diagdata <<"G4InclLightIonInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl;
         
-        diagdata <<"G4InclLightIonInterface: Target A:  " << targetA << G4endl;
-        diagdata <<"G4InclLightIonInterface: Target Z:  " << targetZ << G4endl;
+        diagdata <<"G4InclLightIonInterface: Target A:  " << calincl->targetA() << G4endl;
+        diagdata <<"G4InclLightIonInterface: Target Z:  " << calincl->targetZ() << G4endl;
       }
     }
@@ -175,22 +280,30 @@
       if(bulletType == proton) {
         aParticleDefinition = G4Proton::ProtonDefinition();
-      }
-      if(bulletType == neutron) {
+      } else if(bulletType == neutron) {
         aParticleDefinition = G4Neutron::NeutronDefinition();
-      }
-      if(bulletType == pionPlus) {
+      } else if(bulletType == pionPlus) {
         aParticleDefinition = G4PionPlus::PionPlusDefinition();
-      }
-      if(bulletType == pionZero) {
+      } else if(bulletType == pionZero) {
         aParticleDefinition = G4PionZero::PionZeroDefinition();
-      }
-      if(bulletType == pionMinus) {
+      } else if(bulletType == pionMinus) {
         aParticleDefinition = G4PionMinus::PionMinusDefinition();
-      }
-
-      cascadeParticle = new G4DynamicParticle();
-      cascadeParticle->SetDefinition(aParticleDefinition);
-      cascadeParticle->Set4Momentum(aTrack.Get4Momentum());
-      theResult.AddSecondary(cascadeParticle);
+      } else if(bulletType == deuteron) {
+        aParticleDefinition = G4Deuteron::DeuteronDefinition();
+      } else if(bulletType == triton) {
+        aParticleDefinition = G4Triton::TritonDefinition();
+      } else if(bulletType == he3) {
+        aParticleDefinition = G4He3::He3Definition();
+      } else if(bulletType == he4) {
+        aParticleDefinition = G4Alpha::AlphaDefinition();
+      } else { // Particle was not recognized. Probably an unsupported particle was given as input
+	aParticleDefinition = 0;
+      }
+
+      if(aParticleDefinition != 0) {
+	cascadeParticle = new G4DynamicParticle();
+	cascadeParticle->SetDefinition(aParticleDefinition);
+	cascadeParticle->Set4Momentum(aTrack.Get4Momentum());
+	result.push_back(cascadeParticle);
+      }
     }
 
@@ -199,6 +312,8 @@
     theResult.SetStatusChange(stopAndKill);
     
-    for(particleI = 0; particleI < varntp->ntrack; particleI++) { // Loop through the INCL4+ABLA output.
+    for(particleI = 0; particleI <= varntp->ntrack; particleI++) { // Loop through the INCL4+ABLA output.
       // Get energy/momentum and construct momentum vector in INCL4 coordinates.
+      //      if(varntp->itypcasc[particleI] == -1) continue; // Avoid nucleons that are part of the spectator
+      if(varntp->avv[particleI] == 0 && varntp->zvv[particleI] == 0) continue;
       momx = varntp->plab[particleI]*std::sin(varntp->tetlab[particleI]*CLHEP::pi/180.0)*std::cos(varntp->philab[particleI]*CLHEP::pi/180.0)*MeV;
       momy = varntp->plab[particleI]*std::sin(varntp->tetlab[particleI]*CLHEP::pi/180.0)*std::sin(varntp->philab[particleI]*CLHEP::pi/180.0)*MeV;
@@ -222,4 +337,6 @@
           new G4DynamicParticle(G4Proton::ProtonDefinition(), momDirection, eKin);
         particleIdentified++;
+	baryonNumberBalanceInINCL -= 1;
+	chargeNumberBalanceInINCL -= 1;
       }
 
@@ -228,4 +345,5 @@
           new G4DynamicParticle(G4Neutron::NeutronDefinition(), momDirection, eKin);
         particleIdentified++;
+	baryonNumberBalanceInINCL -= 1;
       }
 
@@ -234,4 +352,5 @@
           new G4DynamicParticle(G4PionPlus::PionPlusDefinition(), momDirection, eKin);
         particleIdentified++;
+	chargeNumberBalanceInINCL -= 1;
       }
 
@@ -240,4 +359,5 @@
           new G4DynamicParticle(G4PionZero::PionZeroDefinition(), momDirection, eKin);
         particleIdentified++;
+	chargeNumberBalanceInINCL -= 0;
       }
 
@@ -246,9 +366,9 @@
           new G4DynamicParticle(G4PionMinus::PionMinusDefinition(), momDirection, eKin);
         particleIdentified++;
+	chargeNumberBalanceInINCL -= -1;
       }
 
       if((varntp->avv[particleI] > 1) && (varntp->zvv[particleI] >= 1)) { // Nucleus fragment
-        G4ParticleDefinition * aIonDef = 0;
-        G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
+	G4ParticleDefinition * aIonDef = 0;
 
         G4int A = G4int(varntp->avv[particleI]);
@@ -273,10 +393,12 @@
 	    new G4DynamicParticle(aIonDef, momDirection, eKin);
 	  particleIdentified++;
+	  baryonNumberBalanceInINCL -= A;
+	  chargeNumberBalanceInINCL -= Z;
 	}
       }
 	
       if(particleIdentified == 1) { // Particle identified properly.
-        cascadeParticle->Set4Momentum(cascadeParticle->Get4Momentum()*=toLabFrame);
-        theResult.AddSecondary(cascadeParticle); // Put data into G4HadFinalState.
+	cascadeParticle->Set4Momentum(cascadeParticle->Get4Momentum()*=toLabFrame);
+        result.push_back(cascadeParticle);
       }
       else { // Particle identification failed.
@@ -292,5 +414,201 @@
     }
 
+    // Spectator nucleus Fermi break-up
+    if(useFermiBreakup && useProjectileSpectator && varntp->masp > 1) {
+      baryonNumberBalanceInINCL -= G4int(varntp->masp);
+      G4double nuclearMass = G4NucleiProperties::GetNuclearMass(G4int(varntp->masp), G4int(varntp->mzsp)) + varntp->exsp * MeV;
+      // Use momentum scaling to compensate for different masses in G4 and INCL:
+      G4double momentumScaling = G4InclUtils::calculate4MomentumScaling(G4int(varntp->masp),
+									G4int(varntp->mzsp),
+									varntp->exsp,
+									varntp->spectatorT,
+									varntp->spectatorP1,
+									varntp->spectatorP2,
+									varntp->spectatorP3);
+      G4LorentzVector p4(momentumScaling * varntp->spectatorP1 * MeV, momentumScaling * varntp->spectatorP2 * MeV,
+			 momentumScaling * varntp->spectatorP3 * MeV,
+			 varntp->spectatorT * MeV + nuclearMass);
+      // Four-momentum, baryon number and charge balance:
+      G4LorentzVector fourMomentumBalance = p4;
+      G4int baryonNumberBalance = G4int(varntp->masp);
+      chargeNumberBalanceInINCL -= G4int(varntp->mzsp);
+      G4int chargeBalance = G4int(varntp->mzsp);
+
+      G4LorentzRotation toFragmentZ;
+      // Assume that Fermi breakup uses Z as the direction of the projectile
+      toFragmentZ.rotateZ(-p4.theta());
+      toFragmentZ.rotateY(-p4.phi());
+      G4LorentzRotation toFragmentLab = toFragmentZ.inverse();
+      p4 *= toFragmentZ;
+      
+      G4LorentzVector p4rest = p4;
+      p4rest.boost(-p4.boostVector());
+      if(verboseLevel > 0) {
+	G4cout <<"Spectator nucleus:" << G4endl;
+	G4cout <<"p4: " << G4endl;
+	G4cout <<" px: " << p4.px() <<" py: " << p4.py() <<" pz: " << p4.pz() << G4endl;
+	G4cout <<" E = " << p4.e() << G4endl;
+	G4cout <<"p4rest: " << G4endl;
+	G4cout <<" px: " << p4rest.px() <<" py: " << p4rest.py() <<" pz: " << p4rest.pz() << G4endl;
+	G4cout <<" E = " << p4rest.e() << G4endl;
+      }
+      G4Fragment theSpectatorNucleus(G4int(varntp->masp), G4int(varntp->mzsp), p4rest);
+      theSpectatorPrecoResult = thePrecoModel->DeExcite(theSpectatorNucleus);
+      if(theSpectatorPrecoResult != 0) {
+      G4ReactionProductVector::iterator fragment;
+      for(fragment = theSpectatorPrecoResult->begin(); fragment != theSpectatorPrecoResult->end(); fragment++) {
+	G4ParticleDefinition *theFragmentDefinition = (*fragment)->GetDefinition();
+
+	if(theFragmentDefinition != 0) {
+	  G4DynamicParticle *theFragment = new G4DynamicParticle(theFragmentDefinition, (*fragment)->GetMomentum());
+	  G4LorentzVector labMomentum = theFragment->Get4Momentum();
+	  labMomentum.boost(p4.boostVector());
+	  labMomentum *= toFragmentLab;
+	  labMomentum *= toLabFrame;
+	  theFragment->Set4Momentum(labMomentum);
+	  fourMomentumBalance -= theFragment->Get4Momentum();
+	  baryonNumberBalance -= theFragmentDefinition->GetAtomicMass();
+	  chargeBalance -= theFragmentDefinition->GetAtomicNumber();
+	  if(verboseLevel > 0) {
+	    G4cout <<"Resulting fragment: " << G4endl;
+	    G4cout <<" kinetic energy = " << theFragment->GetKineticEnergy() / MeV << " MeV" << G4endl;
+	    G4cout <<" momentum = " << theFragment->GetMomentum().mag() / MeV << " MeV" << G4endl;
+	  }
+	  theResult.AddSecondary(theFragment);
+	} else {
+	  G4cout <<"G4InclCascadeInterface: Error. Fragment produced by Fermi break-up does not exist." << G4endl;
+	  G4cout <<"Resulting fragment: " << G4endl;
+	  G4cout <<" momentum = " << (*fragment)->GetMomentum().mag() / MeV << " MeV" << G4endl;
+	}
+      }
+      delete theSpectatorPrecoResult;
+      theSpectatorPrecoResult = 0;
+
+      if(verboseLevel > 1 && std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 
+	G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl;
+	G4cout <<"Magnitude: " << fourMomentumBalance.mag() / MeV << " MeV" << G4endl;
+	G4cout <<"Vector components (px, py, pz, E) = ("
+	       << fourMomentumBalance.px() << ", "
+	       << fourMomentumBalance.py() << ", "
+	       << fourMomentumBalance.pz() << ", "
+	       << fourMomentumBalance.e() << ")" << G4endl;
+      }
+      if(baryonNumberBalance != 0 && verboseLevel > 1) {
+	G4cout <<"Baryon number balance after remnant nucleus Fermi break-up: " << baryonNumberBalance << G4endl;
+      }
+      if(chargeBalance != 0 && verboseLevel > 1) {
+	G4cout <<"Charge balance after remnant nucleus Fermi break-up: " << chargeBalance << G4endl;
+      }
+      }
+    }
+
+    // Finally do Fermi break-up if needed
+    if(varntp->massini > 0) {
+      baryonNumberBalanceInINCL -= G4int(varntp->massini);
+      chargeNumberBalanceInINCL -= G4int(varntp->mzini);
+      // Call Fermi Break-up
+      G4double nuclearMass = G4NucleiProperties::GetNuclearMass(G4int(varntp->massini), G4int(varntp->mzini)) + varntp->exini * MeV;
+      G4LorentzVector fragmentMomentum(varntp->pxrem * MeV, varntp->pyrem * MeV, varntp->pzrem * MeV,
+				       varntp->erecrem * MeV + nuclearMass);
+      G4double momentumScaling = G4InclUtils::calculate4MomentumScaling(G4int(varntp->massini), G4int(varntp->mzini),
+									varntp->exini,
+									varntp->erecrem,
+									varntp->pxrem,
+									varntp->pyrem,
+									varntp->pzrem);
+      G4LorentzVector p4(momentumScaling * varntp->pxrem * MeV, momentumScaling * varntp->pyrem * MeV,
+			 momentumScaling * varntp->pzrem * MeV,
+      			 varntp->erecrem + nuclearMass);
+
+      // For four-momentum, baryon number and charge conservation check:
+      G4LorentzVector fourMomentumBalance = p4;
+      G4int baryonNumberBalance = G4int(varntp->massini);
+      G4int chargeBalance = G4int(varntp->mzini);
+
+      G4LorentzRotation toFragmentZ;
+      toFragmentZ.rotateZ(-p4.theta());
+      toFragmentZ.rotateY(-p4.phi());
+      G4LorentzRotation toFragmentLab = toFragmentZ.inverse();
+      p4 *= toFragmentZ;
+
+      G4LorentzVector p4rest = p4;
+      p4rest.boost(-p4.boostVector());
+      if(verboseLevel > 0) {
+	G4cout <<"Cascade remnant nucleus:" << G4endl;
+	G4cout <<"p4: " << G4endl;
+	G4cout <<" px: " << p4.px() <<" py: " << p4.py() <<" pz: " << p4.pz() << G4endl;
+	G4cout <<" E = " << p4.e() << G4endl;
+
+	G4cout <<"p4rest: " << G4endl;
+	G4cout <<" px: " << p4rest.px() <<" py: " << p4rest.py() <<" pz: " << p4rest.pz() << G4endl;
+	G4cout <<" E = " << p4rest.e() << G4endl;
+      }
+
+      G4Fragment theCascadeRemnant(G4int(varntp->massini), G4int(varntp->mzini), p4rest);
+      thePrecoResult = thePrecoModel->DeExcite(theCascadeRemnant);
+      if(thePrecoResult != 0) {
+      G4ReactionProductVector::iterator fragment;
+      for(fragment = thePrecoResult->begin(); fragment != thePrecoResult->end(); fragment++) {
+	G4ParticleDefinition *theFragmentDefinition = (*fragment)->GetDefinition();
+
+	if(theFragmentDefinition != 0) {
+	  G4DynamicParticle *theFragment = new G4DynamicParticle(theFragmentDefinition, (*fragment)->GetMomentum());
+	  G4LorentzVector labMomentum = theFragment->Get4Momentum();
+	  labMomentum.boost(p4.boostVector());
+	  labMomentum *= toFragmentLab;
+	  labMomentum *= toLabFrame;
+	  theFragment->Set4Momentum(labMomentum);
+	  fourMomentumBalance -= theFragment->Get4Momentum();
+	  baryonNumberBalance -= theFragmentDefinition->GetAtomicMass();
+	  chargeBalance -= theFragmentDefinition->GetAtomicNumber();
+	  if(verboseLevel > 0) {
+	    G4cout <<"Resulting fragment: " << G4endl;
+	    G4cout <<" kinetic energy = " << theFragment->GetKineticEnergy() / MeV << " MeV" << G4endl;
+	    G4cout <<" momentum = " << theFragment->GetMomentum().mag() / MeV << " MeV" << G4endl;
+	  }
+	  theResult.AddSecondary(theFragment);
+	} else {
+	  G4cout <<"G4InclCascadeInterface: Error. Fragment produced by Fermi break-up does not exist." << G4endl;
+	  G4cout <<"Resulting fragment: " << G4endl;
+	  G4cout <<" momentum = " << (*fragment)->GetMomentum().mag() / MeV << " MeV" << G4endl;
+	}
+      }
+      delete thePrecoResult;
+      thePrecoResult = 0;
+
+      if(verboseLevel > 1 && std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 
+	G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl;
+	G4cout <<"Magnitude: " << fourMomentumBalance.mag() / MeV << " MeV" << G4endl;
+	G4cout <<"Vector components (px, py, pz, E) = ("
+	       << fourMomentumBalance.px() << ", "
+	       << fourMomentumBalance.py() << ", "
+	       << fourMomentumBalance.pz() << ", "
+	       << fourMomentumBalance.e() << ")" << G4endl;
+      }
+      if(baryonNumberBalance != 0 && verboseLevel > 1) {
+	G4cout <<"Baryon number balance after remnant nucleus Fermi break-up: " << baryonNumberBalance << G4endl;
+      }
+      if(chargeBalance != 0 && verboseLevel > 1) {
+	G4cout <<"Charge balance after remnant nucleus Fermi break-up: " << chargeBalance << G4endl;
+      }
+      }
+    }
+
     varntp->ntrack = 0; // Clean up the number of generated particles in the event.
+
+    if(baryonNumberBalanceInINCL != 0 && verboseLevel > 1) {
+      G4cout <<"Event " << eventNumber <<": G4InclLightIonInterface: Baryon number conservation problem in INCL detected!" << G4endl;
+      G4cout <<"Baryon number balance: " << baryonNumberBalanceInINCL << G4endl;
+      if(baryonNumberBalanceInINCL < 0) {
+	G4cout <<"Event " << eventNumber <<": Too many outcoming baryons!" << G4endl;
+      } else if(baryonNumberBalanceInINCL > 0) {
+	G4cout <<"Event " << eventNumber <<": Too few outcoming baryons!" << G4endl;
+      }
+    }
+
+    if(chargeNumberBalanceInINCL != 0 && verboseLevel > 1) {
+      G4cout <<"Event " << eventNumber <<": G4InclLightIonInterface: Charge number conservation problem in INCL detected!" << G4endl;
+      G4cout <<"Event " << eventNumber <<": Charge number balance: " << chargeNumberBalanceInINCL << G4endl;
+    }
   }
   /**
@@ -304,5 +622,5 @@
     cascadeParticle = new G4DynamicParticle(theTableOfParticles->FindParticle(aTrack.GetDefinition()), aTrack.Get4Momentum());
 
-    theResult.AddSecondary(cascadeParticle);
+    result.push_back(cascadeParticle);
 
     if(verboseLevel > 1) {
@@ -317,5 +635,5 @@
 	G4cout <<"G4InclLightIonInterface: Unknown bullet type" << G4endl;      
 	G4cout <<"Bullet particle name: " << cascadeParticle->GetDefinition()->GetParticleName() << G4endl;
-      }      
+      }
       if(verboseLevel > 3) {
         diagdata <<"G4InclLightIonInterface: Unknown bullet type" << G4endl;      
@@ -324,26 +642,27 @@
     }
 
-    if((targetA == 1) && (targetZ == 1)) { // Unsupported target
+    if((calincl->targetA() == 1) && (calincl->targetZ() == 1)) { // Unsupported target
       if(verboseLevel > 1) {
 	G4cout <<"Unsupported target: " << G4endl;
-	G4cout <<"Target A: " << targetA << G4endl;
-	G4cout <<"TargetZ: " << targetZ << G4endl;
+	G4cout <<"Target A: " << calincl->targetA() << G4endl;
+	G4cout <<"TargetZ: " << calincl->targetZ() << G4endl;
       }
       if(verboseLevel > 3) {
         diagdata <<"Unsupported target: " << G4endl;
-        diagdata <<"Target A: " << targetA << G4endl;
-       diagdata <<"TargetZ: " << targetZ << G4endl;
-      }
-    }
-
-    if(bulletE < 100) { // INCL does not support E < 100 MeV.
+        diagdata <<"Target A: " << calincl->targetA() << G4endl;
+	diagdata <<"TargetZ: " << calincl->targetZ() << G4endl;
+      }
+    }
+
+    if(calincl->bulletE() < 100) { // INCL does not support E < 100 MeV.
       if(verboseLevel > 1) {
-	G4cout <<"Unsupported bullet energy: " << bulletE << " MeV. (Lower limit is 100 MeV)." << G4endl;
+	G4cout <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;
 	G4cout <<"WARNING: Returning the original bullet with original energy back to Geant4." << G4endl;
       }
       if(verboseLevel > 3) {
-        diagdata <<"Unsupported bullet energy: " << bulletE << " MeV. (Lower limit is 100 MeV)." << G4endl;
-      }
-    }
+        diagdata <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;
+      }
+    }
+
     if(verboseLevel > 3) {
       diagdata <<"WARNING: returning the original bullet with original energy back to Geant4." << G4endl;
@@ -351,4 +670,22 @@
   }
 
+  // Finally copy the accumulated secondaries into the result collection:
+  G4ThreeVector boostVector = aTrack.Get4Momentum().boostVector();
+  G4LorentzRotation boostBack = toBreit.inverse();
+
+  for(std::vector<G4DynamicParticle*>::iterator i = result.begin(); i != result.end(); ++i) {
+    // If the calculation was performed in inverse kinematics we have to
+    // convert the result back...
+    if(calincl->isInverseKinematics()) {
+      G4LorentzVector mom = (*i)->Get4Momentum();
+      mom.setPz(-1.0 * mom.pz()); // Reverse the z-component of the momentum vector
+      mom *= boostBack;
+      (*i)->Set4Momentum(mom);
+    }
+    theResult.AddSecondary((*i));
+  }
+
+  delete calincl;
+  calincl = 0;
   return &theResult;
 } 
