Changeset 1347 for trunk/source/processes/hadronic/models/incl

Timestamp:

Dec 22, 2010, 3:52:27 PM (14 years ago)

Author:

garnier

Message:

geant4 tag 9.4

Location:

trunk/source/processes/hadronic/models/incl

Files:

• GNUmakefile (modified) (2 diffs)
• History (modified) (2 diffs)
• include/G4Abla.hh (modified) (4 diffs)
• include/G4AblaDataDefs.hh (modified) (5 diffs)
• include/G4Incl.hh (modified) (2 diffs)
• include/G4InclAblaCascadeInterface.hh (modified) (2 diffs)
• include/G4InclAblaDataFile.hh (modified) (2 diffs)
• include/G4InclCascadeInterface.hh (modified) (7 diffs)
• include/G4InclDataDefs.hh (modified) (31 diffs)
• include/G4InclLightIonInterface.hh (modified) (6 diffs)
• include/G4InclRandomNumbers.hh (modified) (2 diffs)
• src/G4Abla.cc (modified) (9 diffs)
• src/G4AblaFission.cc (modified) (3 diffs)
• src/G4AblaFissionSimfis18.cc (modified) (3 diffs)
• src/G4Incl.cc (modified) (24 diffs)
• src/G4InclAblaCascadeInterface.cc (modified) (3 diffs)
• src/G4InclAblaDataFile.cc (modified) (3 diffs)
• src/G4InclAblaLightIonInterface.cc (modified) (6 diffs)
• src/G4InclCascadeInterface.cc (modified) (22 diffs)
• src/G4InclLightIonInterface.cc (modified) (22 diffs)

trunk/source/processes/hadronic/models/incl/GNUmakefile

@@ -1 @@
-# $Id: GNUmakefile,v 1.11 2010/10/26 02:47:58 kaitanie Exp $
+# $Id: GNUmakefile,v 1.13 2010/11/13 00:08:36 kaitanie Exp $
 # -----------------------------------------------------------
 # GNUmakefile for hadronic library.  Gabriele Cosmo, 18/9/96.
@@ -36 +36 @@
             -I$(G4BASE)/processes/hadronic/models/management/include \
             -I$(G4BASE)/processes/hadronic/models/util/include \
+            -I$(G4BASE)/processes/hadronic/models/de_excitation/util/include \
             -I$(G4BASE)/processes/hadronic/models/de_excitation/evaporation/include \
             -I$(G4BASE)/processes/hadronic/models/de_excitation/fermi_breakup/include \
+            -I$(G4BASE)/processes/hadronic/models/de_excitation/photon_evaporation/include \
+            -I$(G4BASE)/processes/hadronic/models/de_excitation/multifragmentation/include \
+            -I$(G4BASE)/processes/hadronic/models/de_excitation/handler/include \
+            -I$(G4BASE)/processes/hadronic/models/de_excitation/management/include \
+            -I$(G4BASE)/processes/hadronic/models/pre_equilibrium/exciton_model/include/ \
             -I$(G4BASE)/particles/management/include \
             -I$(G4BASE)/particles/leptons/include \

trunk/source/processes/hadronic/models/incl/History

@@ -4 +4 @@
      Geant4 - an Object-Oriented Toolkit for Physics Simulation
      ==========================================================
-$Id: History,v 1.32 2010/10/29 06:54:23 gunter Exp $
+$Id: History,v 1.37 2010/11/17 20:19:09 kaitanie Exp $
 ---------------------------------------------------------------------
 
@@ -17 +17 @@
    ---------------------------------------------------------------
 
+17 November 2010 - Pekka Kaitaniemi (hadr-incl-V09-03-10)
+---------------------------------------------------------
+- Fixed several issues reported by Coverity:
+  o Fix: Fragment vector and Fermi break-up related emory leaks in
+  INCL/ABLA interfaces
+  o Initialize INCL internal variables properly
+  o Check array boundaries in datafile reader
+
+16 November 2010 - Gabriele Cosmo (hadr-incl-V09-03-09)
+-------------------------------------------------------
+- More minor fixes from Coverity reports...
+
+12 November 2010 - Pekka Kaitaniemi (hadr-incl-V09-03-08)
+---------------------------------------------------------
+- Fixed several minor variable initialization issues reported by Coverity
+
+11 November 2010 - Pekka Kaitaniemi (hadr-incl-V09-03-07)
+---------------------------------------------------------
+- Fixed some minor variable initialization issues reported by Coverity
+- Updated interfaces:
+  o INCL + built-in ABLA de-excitation: G4InclAblaCascadeInterface and
+  G4InclAblaLightIonInterface
+  o INCL + PreCompound model: G4InclCascadeInterface and
+  G4InclLightIonInterface
+
+03 November 2010 - Pekka Kaitaniemi (hadr-incl-V09-03-06)
+---------------------------------------------------------
+- Fixed insufficient array index safeguard in ABLA
+- Silenced debugging output printed by G4InclAblaCascadeInterface
+
 29 October 2010 - Gunter Folger (hadr-incl-V09-03-05)
+-----------------------------------------------------
 - Fixed several more compilation warnings for conversion of double to int
    in  G4Incl.cc, G4InclAblaCascadeInterface.cc, and G4InclAblaLightIonInterface.cc.

trunk/source/processes/hadronic/models/incl/include/G4Abla.hh

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4Abla.hh,v 1.13 2010/10/26 02:47:59 kaitanie Exp $ 
+// $Id: G4Abla.hh,v 1.14 2010/11/17 20:19:09 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -50 +50 @@
 public:
   /**
-   * Basic constructor.
-   */
-  G4Abla();
-
-  /**
    * This constructor is used by standalone test driver and the Geant4 interface.
    *
@@ -62 +57 @@
    */
   G4Abla(G4Hazard *aHazard, G4Volant *aVolant, G4VarNtp *aVarntp);
-
-  /**
-   * Constructor that is to be used only for testing purposes.
-   * @param aHazard random seeds
-   * @param aVolant data structure for ABLA output   
-   */
-  G4Abla(G4Hazard *hazard, G4Volant *volant);
 
   /**
@@ -167 +155 @@
               G4double *probp_par, G4double *probn_par, G4double *proba_par, 
               G4double *probf_par, G4double *ptotl_par, G4double *sn_par, 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 *ecp_par,G4double *eca_par, G4double *bp_par, G4double *ba_par, G4int, G4int inum, G4int itest);
 
   /**

trunk/source/processes/hadronic/models/incl/include/G4AblaDataDefs.hh

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4AblaDataDefs.hh,v 1.11 2010/10/26 02:47:59 kaitanie Exp $ 
+// $Id: G4AblaDataDefs.hh,v 1.13 2010/11/13 00:08:36 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -95 +95 @@
    * 
    */
-  G4Ald() {};
+  G4Ald()
+    :av(0.0), as(0.0), ak(0.0), optafan(0.0)
+  {};
   ~G4Ald() {};
   
@@ -141 +143 @@
 
 public:
-  G4Fiss() {};
+  G4Fiss()
+    :akap(0.0), bet(0.0), homega(0.0), koeff(0.0), ifis(0.0),
+     optshp(0), optxfis(0), optles(0), optcol(0)
+  {};
   ~G4Fiss() {};
   
@@ -171 +176 @@
 
 public:
-  G4Opt() {};
+  G4Opt()
+    :optemd(0), optcha(0), eefac(0.0)
+  {};
   ~G4Opt() {};
   
@@ -182 +189 @@
 class G4Eenuc {
 public:
-  G4Eenuc() {};
+  G4Eenuc() {
+    for(G4int i = 0; i < EENUCSIZE; ++i) {
+      she[i] = 0.0;
+    }
+    for(G4int i = 0; i < XHESIZE; ++i) {
+      for(G4int j = 0; j < EENUCSIZE; ++j) {
+        xhe[i][j] = 0.0;
+      }
+    }
+  };
   ~G4Eenuc() {};
   

trunk/source/processes/hadronic/models/incl/include/G4Incl.hh

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4Incl.hh,v 1.18 2010/10/26 02:47:59 kaitanie Exp $ 
+// $Id: G4Incl.hh,v 1.19 2010/11/13 00:08:36 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -688 +688 @@
    private:
 
+  /*
+   * (Re)Initialize INCL internal variables
+   */
+  void clearState();
+
   /**
    * Random seeds for INCL4 internal random number generators.

trunk/source/processes/hadronic/models/incl/include/G4InclAblaCascadeInterface.hh

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4InclAblaCascadeInterface.hh,v 1.9 2010/10/26 02:47:59 kaitanie Exp $ 
+// $Id: G4InclAblaCascadeInterface.hh,v 1.11 2010/11/13 00:08:36 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -91 +91 @@
  *
  * @see G4InclAblaLightIonInterface
+ * @see G4InclCascadeInterface
  */
 

trunk/source/processes/hadronic/models/incl/include/G4InclAblaDataFile.hh

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4InclAblaDataFile.hh,v 1.3 2010/06/14 16:10:01 gcosmo Exp $ 
+// $Id: G4InclAblaDataFile.hh,v 1.4 2010/11/17 20:19:09 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -52 +52 @@
 private:
   G4int verboseLevel;
-  
-  G4String *dataPath;
 };
 

trunk/source/processes/hadronic/models/incl/include/G4InclCascadeInterface.hh

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4InclCascadeInterface.hh,v 1.6 2010/10/26 02:47:59 kaitanie Exp $ 
+// $Id: G4InclCascadeInterface.hh,v 1.8 2010/11/13 00:08:36 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -59 +59 @@
 #include "G4AblaDataDefs.hh"
 #include "G4Incl.hh"
-#include "G4InclInput.hh"
+
+// Geant4 de-excitation
+#include "G4ExcitationHandler.hh"
+#include "G4PreCompoundModel.hh"
 
 #include <fstream>
@@ -67 +70 @@
 
 /**
+ * <h1>INCL intra-nuclear cascade with Geant4 PreCompound for de-excitation</h1>
+ *
  * Interface for INCL. This interface handles basic hadron
  * bullet particles (protons, neutrons, pions).
- * @see G4InclAblaLightIonInterface
+ *
+ * Example usage in case of protons:
+ * @code
+ * G4InclCascadeInterface* inclModel = new G4InclCascadeInterface;
+ * inclModel -> SetMinEnergy(0.0 * MeV); // Set the energy limits
+ * inclModel -> SetMaxEnergy(3.0 * GeV);
+ *
+ * G4ProtonInelasticProcess* protonInelasticProcess = new G4ProtonInelasticProcess(); 
+ * G4ProtonInelasticCrossSection* protonInelasticCrossSection =  new G4ProtonInelasticCrossSection(); 
+ *
+ * protonInelasticProcess -> RegisterMe(inclModel);
+ * protonInelasticProcess -> AddDataSet(protonInelasticCrossSection);
+ *
+ * particle = G4Proton::Proton();
+ * processManager = particle -> GetProcessManager();
+ * processManager -> AddDiscreteProcess(protonInelasticProcess);
+ * @endcode
+ * The same setup procedure is needed for neutron and pion inelastic processes
+ * as well.
+ *
+ * @see G4InclLightIonInterface
  */
 
@@ -78 +103 @@
    * Basic constructor.
    */
-  G4InclCascadeInterface();
+  G4InclCascadeInterface(const G4String& name = "INCL Cascade with Geant4 PreCompound");
 
   
@@ -94 +119 @@
 
   /**
-   * Main method to apply the INCL/ABLA physics model.
+   * Main method to apply the INCL physics model.
    * @param aTrack the projectile particle
    * @param theNucleus target nucleus
-   * @return the output of the INCL/ABLA physics model
+   * @return the output of the INCL physics model
    */
   G4HadFinalState* ApplyYourself(const G4HadProjectile& aTrack,  G4Nucleus& theNucleus); 
@@ -108 +133 @@
   G4Hazard *hazard; // The random seeds used by INCL.
   G4VarNtp *varntp;
-  G4InclInput *inclInput;
+  G4InclInput *calincl;
   G4Ws *ws;
   G4Mat *mat;
   G4Incl *incl;
-  
+
   G4HadFinalState theResult;  
   ofstream diagdata;
@@ -119 +144 @@
   G4double previousTargetA;
   G4double previousTargetZ;
+
+  G4ExcitationHandler *theExcitationHandler;
+  G4PreCompoundModel *thePrecoModel;
 };
 
-#endif // G4INCLCASCADEINTERFACE_H
+#endif // G4INCLABLACASCADEINTERFACE_H

trunk/source/processes/hadronic/models/incl/include/G4InclDataDefs.hh

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4InclDataDefs.hh,v 1.10 2010/10/28 15:35:50 gcosmo Exp $ 
+// $Id: G4InclDataDefs.hh,v 1.12 2010/11/17 20:19:09 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -63 +63 @@
 class G4QuadvectProjo {
 public:
-  G4QuadvectProjo() {};
+  G4QuadvectProjo() {
+    for(G4int i = 0; i < max_a_proj; ++i) {
+      eps_c[i] = 0.0;
+      t_c[i] = 0.0;
+      p3_c[i] = 0.0;
+      p1_s[i] = 0.0;
+      p2_s[i] = 0.0;
+      p3_s[i] = 0.0;
+    }
+  };
+
   ~G4QuadvectProjo() {};
 
@@ -73 +83 @@
 class G4VBe {
 public:
-  G4VBe() {};
+  G4VBe()
+    :ia_be(0), iz_be(0),
+     rms_be(0.0), pms_be(0.0), bind_be(0.0)
+  { };
+
   ~G4VBe() {};
 
@@ -87 +101 @@
   G4InclProjSpect() {
     //    G4cout <<"Projectile spectator data structure created!" << G4endl;
+    clear();
   };
   ~G4InclProjSpect() {};
@@ -95 +110 @@
     a_projspec = 0;
     z_projspec = 0;
-    m_projspec = 0.0;
+    t_projspec = 0.0;
     ex_projspec = 0.0;
     p1_projspec = 0.0;
     p2_projspec = 0.0;
     p3_projspec = 0.0;
+    m_projspec = 0.0;
   };
 
@@ -108 +124 @@
 };
 
-#define FSIZE 15
-/**
- * Initial values of a hadronic cascade problem.
- */
-class G4Calincl {
-public:
-  G4Calincl() {
-    isExtendedProjectile = false;
-  };
-
-  G4Calincl(const G4HadProjectile &aTrack, const G4Nucleus &theNucleus, G4bool inverseKinematics = false) {
-    for(int i = 0; i < 15; i++) {
-      f[i] = 0.0; // Initialize INCL input data
-    }
-
-    usingInverseKinematics = inverseKinematics;
-
-    f[0] = theNucleus.GetA_asInt(); // Target mass number
-    f[1] = theNucleus.GetZ_asInt(); // Target charge number
-    f[6] = getBulletType(aTrack.GetDefinition()); // Projectile type (INCL particle code)
-    f[2] = aTrack.GetKineticEnergy() / MeV; // Projectile energy (total, in MeV)
-    f[5] = 1.0; // Time scaling
-    f[4] = 45.0; // Nuclear potential
-    setExtendedProjectileInfo(aTrack.GetDefinition());
-    //    G4cout <<"Projectile type = " << f[6] << G4endl;
-    //    G4cout <<"Energy = " << f[2] << G4endl;
-    //    G4cout <<"Target A = " << f[0] << " Z = " << f[1] << G4endl;
-  };
-
-  ~G4Calincl() {};
-  
-  static void printProjectileTargetInfo(const G4HadProjectile &aTrack, const G4Nucleus &theNucleus) {
-    G4cout <<"Projectile = " << aTrack.GetDefinition()->GetParticleName() << G4endl;
-    G4cout <<"    four-momentum: " << aTrack.Get4Momentum() << G4endl;
-    G4cout <<"Energy = " << aTrack.GetKineticEnergy() / MeV << G4endl;
-    G4cout <<"Target A = " << theNucleus.GetA_asInt() << " Z = " << theNucleus.GetZ_asInt() << G4endl;
-  }
-
-  static G4bool canUseInverseKinematics(const G4HadProjectile &aTrack, const G4Nucleus &theNucleus) {
-    G4int targetA = theNucleus.GetA_asInt();
-    const G4ParticleDefinition *projectileDef = aTrack.GetDefinition();
-    G4int projectileA = projectileDef->GetAtomicMass();
-    //    G4int projectileZ = projectileDef->GetAtomicNumber();
-    if(targetA > 0 && targetA < 18 && (projectileDef != G4Proton::Proton() &&
-                                       projectileDef != G4Neutron::Neutron() &&
-                                       projectileDef != G4PionPlus::PionPlus() &&
-                                       projectileDef != G4PionZero::PionZero() &&
-                                       projectileDef != G4PionMinus::PionMinus()) &&
-       projectileA > 1) {
-      return true;
-    } else {
-      return false;
-    }
-  }
-
-  G4double bulletE() {
-    return f[2];
-  }
-
-  G4int getBulletType() {
-    return G4int(f[6]);
-  };
-
-  void setExtendedProjectileInfo(const G4ParticleDefinition *pd) {
-    if(getBulletType(pd) == -666) {
-      extendedProjectileA = pd->GetAtomicMass();
-      extendedProjectileZ = pd->GetAtomicNumber();
-      isExtendedProjectile = true;
-    } else {
-      isExtendedProjectile = false;
-    }
-  };
-
-  G4int getBulletType(const G4ParticleDefinition *pd) {
-    G4ParticleTable *pt = G4ParticleTable::GetParticleTable();
-
-    if(pd == G4Proton::Proton()) {
-      return 1;
-    } else if(pd == G4Neutron::Neutron()) {
-      return 2;
-    } else if(pd == G4PionPlus::PionPlus()) {
-      return 3;
-    } else if(pd == G4PionMinus::PionMinus()) {
-      return 5;
-    } else if(pd == G4PionZero::PionZero()) {
-      return 4;
-    } else if(pd == G4Deuteron::Deuteron()) {
-      return 6;
-    } else if(pd == G4Triton::Triton()) {
-      return 7;
-    } else if(pd == G4He3::He3()) {
-      return 8;
-    } else if(pd == G4Alpha::Alpha()) {
-      return 9;
-    } else if(pd == pt->GetIon(6, 12, 0.0)) { // C12 special case. This should be phased-out in favor of "extended projectile"
-      return -12;
-    } else { // Is this extended projectile?
-      G4int A = pd->GetAtomicMass();
-      G4int Z = pd->GetAtomicNumber();
-      if(A > 4 && A <= 16 && Z > 2 && Z <= 8) { // Ions from Lithium to Oxygen
-        return -666; // Code of an extended projectile
-      }
-    }
-    G4cout <<"Error! Projectile " << pd->GetParticleName() << " not defined!" << G4endl;
-    return 0;
-  };
-
-  G4bool isInverseKinematics() {
-    return usingInverseKinematics;
-  };
-
-  G4double targetA() { return f[0]; };
-  G4double targetZ() { return f[1]; };
-
-  G4int extendedProjectileA;
-  G4int extendedProjectileZ;
-  G4bool isExtendedProjectile;
-
-  /**
-   * Here f is an array containing the following initial values:
-   * - f[0] : target mass number
-   * - f[1] : target charge number
-   * - f[2] : bullet energy
-   * - f[3] : minimum proton energy to leave the target (default: 0.0)
-   * - f[4] : nuclear potential (default: 45.0 MeV)
-   * - f[5] : time scale (default: 1.0)
-   * - f[6] : bullet type (1: proton, 2: neutron, 3: pi+, 4: pi0 5: pi-, 6:H2, 7: H3, 8: He3, 9: He4
-   * - f[7] : minimum neutron energy to leave the target (default: 0.0)
-   * - f[8] : target material identifier (G4Mat)
-   * - f[9] : not used
-   * - f[10] : not used
-   * - f[11] : not used
-   * - f[12] : not used
-   * - f[13] : not used
-   * - f[14] : not used
-   */
-  G4double f[FSIZE];
-
-  /**
-   * Number of events to be processed.
-   */
-  G4int icoup;
-
-  G4bool usingInverseKinematics;
-};
-
 #define IGRAINESIZE 19
 /**
@@ -262 +132 @@
 class G4Hazard{
 public:
-  G4Hazard() {};
+  G4Hazard() {
+    ial = 0;
+    for(G4int i = 0; i < IGRAINESIZE; ++i) {
+      igraine[i] = 0;
+    }
+  };
+
   ~G4Hazard() {};
 
@@ -283 +159 @@
 class G4Mat {
 public:
-  G4Mat() { };
+  G4Mat() {
+    nbmat = 0;
+
+    for(G4int i = 0; i < MATSIZE; ++i) {
+      zmat[i] = 0;
+      amat[i] = 0;
+    }
+
+    for(G4int i = 0; i < MATGEOSIZE; ++i) {
+      for(G4int j = 0; j < MATSIZE; ++j) {
+        bmax_geo[i][j] = 0;
+      }
+    }
+};
+
   ~G4Mat() { };
 
@@ -313 +203 @@
 class G4LightGausNuc {
 public:
-  G4LightGausNuc() {};
+  G4LightGausNuc() {
+    for(G4int i = 0; i < LGNSIZE; ++i) {
+      rms1t[i] = 0.0;
+      pf1t[i] = 0.0;
+      pfln[i] = 0.0;
+      tfln[i] = 0.0;
+      vnuc[i] = 0.0;
+    }
+  };
+
   ~G4LightGausNuc() {};
   
@@ -329 +228 @@
 class G4LightNuc {
 public:
-  G4LightNuc() {};
+  G4LightNuc() {
+    for(G4int i = 0; i < LNSIZE; ++i) {
+      r[i] = 0.0;
+      a[i] = 0.0;
+    }
+  };
+
   ~G4LightNuc() {};
 
@@ -350 +255 @@
 class G4Saxw {
 public:
-  G4Saxw() {};
+  G4Saxw() {
+    for(G4int i = 0; i < SAXWROWS; ++i) {
+      for(G4int j = 0; j < SAXWCOLS; ++j) {
+        x[i][j] = 0.0;
+        y[i][j] = 0.0;
+        s[i][j] = 0.0;
+      }
+    }
+    imat = 0; n = 0; k = 0;
+  };
+
   ~G4Saxw() {};
   
@@ -391 +306 @@
   G4Ws() {
     fneck = 0.0;
-  };
+    r0 = 0.0;
+    adif = 0.0;
+    rmaxws = 0.0;
+    drws = 0.0;
+    nosurf = 0.0;
+    xfoisa = 0.0;
+    bmax = 0.0;
+    npaulstr = 0.0;
+  };
+
   ~G4Ws() {};
   
@@ -453 +377 @@
 class G4Dton {
 public:
-  G4Dton() {};
+  G4Dton() {
+    fn = 0.0;
+    for(G4int i = 0; i < DTONSIZE; ++i) {
+      c[i] = 0.0;
+      d[i] = 0.0;
+    }
+  };
+
   ~G4Dton() {};
   
@@ -469 +400 @@
 class G4Spl2 {
 public:
-  G4Spl2() {};
+  G4Spl2() {
+    for(G4int i = 0; i < SPL2SIZE; ++i) {
+      x[i] = 0.0; y[i] = 0.0;
+      a[i] = 0.0; b[i] = 0.0; c[i] = 0.0;
+    }
+    n = 0;
+  };
+
   ~G4Spl2() {};
   
@@ -491 +429 @@
 class G4Bl1 {
 public:
-  G4Bl1() {};
+  G4Bl1() {
+    ta = 0.0;
+    for(G4int i = 0; i < BL1SIZE; ++i) {
+      p1[i] = 0.0; p2[i] = 0.0; p3[i] = 0.0; eps[i] = 0.0;
+      ind1[i] = 0; ind2[i] = 0;
+    }
+  };
+
   ~G4Bl1() {};
   
@@ -510 +455 @@
 };
 
-#define BL2CROISSIZE 19900
-#define BL2INDSIZE 19900
+#define BL2SIZE 19900
 /**
  * 
@@ -517 +461 @@
 class G4Bl2 {
 public:
-  G4Bl2() {};
+  G4Bl2() {
+    k = 0;
+    for(G4int i = 0; i < BL2SIZE; ++i) {
+      crois[i] = 0.0;
+      ind[i] = 0;
+      jnd[i] = 0;
+    }    
+  };
+
   ~G4Bl2() {};
   
@@ -533 +485 @@
    * 
    */
-  G4double crois[BL2CROISSIZE];
+  G4double crois[BL2SIZE];
 
   /**
@@ -543 +495 @@
    *
    */
-  G4int ind[BL2INDSIZE];
-
-  /**
-   *
-   */
-  G4int jnd[BL2INDSIZE];
+  G4int ind[BL2SIZE];
+
+  /**
+   *
+   */
+  G4int jnd[BL2SIZE];
 };
 
@@ -558 +510 @@
 class G4Bl3 {
 public:
-  G4Bl3() {};
+  G4Bl3() {
+    r1 = 0.0; r2 = 0.0;
+    ia1 = 0; ia2 = 0;
+    rab2 = 0.0;
+
+    for(G4int i = 0; i < BL3SIZE; ++i) {
+      x1[i] = 0.0; x2[i] = 0.0; x3[i] = 0.0;
+    }
+  };
+
   ~G4Bl3() {};
   
@@ -599 +560 @@
 class G4Bl4 {
 public:
-  G4Bl4() {};
+  G4Bl4()
+    :tmax5(0.0)
+  {};
+
   ~G4Bl4() {};
 
@@ -615 +579 @@
 class G4Bl5 {
 public:
-  G4Bl5() {};
+  G4Bl5() {
+    for(G4int i = 0; i < BL5SIZE; ++i) {
+      tlg[i] = 0.0;
+      nesc[i] = 0;
+    }
+  };
+
   ~G4Bl5() {};
   
@@ -634 +604 @@
 class G4Bl6 {
 public:
-  G4Bl6() {};
+  G4Bl6()
+    :xx10(0.0), isa(0.0)
+  {};
+
   ~G4Bl6() {};
   
@@ -653 +626 @@
 class G4Bl8 {
 public:
-  G4Bl8() {};
+  G4Bl8()
+    :rathr(0.0), ramass(0.0)
+  {};
+
   ~G4Bl8() {};
 
@@ -677 +653 @@
     l1 = 0;
     l2 = 0;
+    for(G4int i = 0; i < BL9SIZE; ++i) {
+      hel[i] = 0.0;
+    }
   };
   ~G4Bl9() {};
@@ -696 +675 @@
 class G4Bl10 {
 public:
-  G4Bl10() {};
+  G4Bl10()
+    :ri4(0.0), rs4(0.0), r2i(0.0), r2s(0.0), pdummy(0.0), pf(0.0)
+  {};
   ~G4Bl10() {};
 
@@ -710 +691 @@
 class G4Kind {
 public:
-  G4Kind() {};
+  G4Kind()
+    :kindf7(0)
+  {};
   ~G4Kind() {};
 
@@ -770 +753 @@
 class G4VarAvat {
 public:
-  G4VarAvat() {};
+  G4VarAvat() {
+    kveux = 0;
+    bavat = 0.0;
+    nopartavat = 0; ncolavat = 0;
+    nb_avat = 0;
+
+    for(G4int i = 0; i < VARSIZE; ++i) {
+      r1_in[i] = 0.0;
+      r1_first_avat[i] = 0.0;
+    }
+
+    for(G4int i = 0; i < VAEPSSIZE; ++i) {
+      epsd[i] = 0.0;
+      eps2[i] = 0.0;
+      eps4[i] = 0.0;
+      eps6[i] = 0.0;
+      epsf[i] = 0.0;
+    }
+
+    for(G4int i = 0; i < VAAVM; ++i) {
+      timeavat[i] = 0.0;
+      l1avat[i] = 0.0;
+      l2avat[i] = 0.0;
+      jpartl1[i] = 0.0;
+      jpartl2[i] = 0.0;
+      del1avat[i] = 0.0;
+      del2avat[i] = 0.0;
+      energyavat[i] = 0.0;
+      bloc_paul[i] = 0.0;
+      bloc_cdpp[i] = 0.0;
+      go_out[i] = 0.0;
+    }
+  };
+
   ~G4VarAvat() {};
 
@@ -822 +838 @@
 class G4VarNtp {
 public:
-  G4VarNtp() {};
+  G4VarNtp() {
+    clear();
+  };
+
   ~G4VarNtp() {};
 
@@ -834 +853 @@
     targetA = 0;
     targetZ = 0;
-    masp = 0.0; mzsp = 0.0; exsp = 0.0;
+    masp = 0.0; mzsp = 0.0; exsp = 0.0; mrem = 0.0;
     // To be deleted?
     spectatorA = 0;
@@ -1190 +1209 @@
 class G4Paul {
 public:
-  G4Paul() {};
+  G4Paul()
+    :ct0(0.0), ct1(0.0), ct2(0.0), ct3(0.0), ct4(0.0), ct5(0.0), ct6(0.0),
+     pr(0.0), pr2(0.0), xrr(0.0), xrr2(0.0),
+     cp0(0.0), cp1(0.0), cp2(0.0), cp3(0.0), cp4(0.0), cp5(0.0), cp6(0.0)
+  {};
+
   ~G4Paul() {};
   

trunk/source/processes/hadronic/models/incl/include/G4InclLightIonInterface.hh

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4InclLightIonInterface.hh,v 1.6 2010/10/26 02:47:59 kaitanie Exp $ 
+// $Id: G4InclLightIonInterface.hh,v 1.8 2010/11/13 00:08:36 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -46 +46 @@
 #define G4INCLLIGHTIONINTERFACE_H 1
 
-#include "G4Nucleon.hh"
-#include "G4Nucleus.hh"
-#include "G4HadronicInteraction.hh"
 #include "G4VIntraNuclearTransportModel.hh"
 #include "G4KineticTrackVector.hh"
@@ -61 +58 @@
 #include "G4Incl.hh"
 
+// Geant4 de-excitation
+#include "G4ExcitationHandler.hh"
+#include "G4PreCompoundModel.hh"
+
 #include <fstream>
 #include <iostream>
@@ -67 +68 @@
 
 /**
- * Interface for INCL. This interface handles basic light ion
- * bullet particles (deuterons, tritons, he3 and alphas).
- * @see G4InclAblaLightIonInterface
+ * Interface for INCL with Geant4 PreCompound de-excitation. This
+ * interface handles basic light ion bullet particles from deuterons
+ * up to carbon-12. @see G4InclAblaLightIonInterface
  */
 
@@ -94 +95 @@
 
   /**
-   * Main method to apply the INCL/ABLA physics model.
+   * Main method to apply the INCL physics model.
    * @param aTrack the projectile particle
    * @param theNucleus target nucleus
-   * @return the output of the INCL/ABLA physics model
+   * @return the output of the INCL physics model
    */
   G4HadFinalState* ApplyYourself(const G4HadProjectile& aTrack,  G4Nucleus& theNucleus); 
@@ -119 +120 @@
   G4double previousTargetA;
   G4double previousTargetZ;
+  G4bool useProjectileSpectator;
+  G4bool useFermiBreakup;
+
+  G4ExcitationHandler *theExcitationHandler;
+  G4PreCompoundModel *thePrecoModel;
 };
 
-#endif // G4INCLLIGHTIONINTERFACE_H
+#endif // G4INCLABLALIGHTIONINTERFACE_H

trunk/source/processes/hadronic/models/incl/include/G4InclRandomNumbers.hh

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4InclRandomNumbers.hh,v 1.6 2010/10/26 02:47:59 kaitanie Exp $
+// $Id: G4InclRandomNumbers.hh,v 1.8 2010/11/13 00:08:36 kaitanie Exp $
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -43 +43 @@
 
 public:
-  G4InclRandomInterface() { }
+  G4InclRandomInterface() {
+    this->seed = 1337; // Default seed, this is never actually used.
+  }
   G4InclRandomInterface(G4long seed) {
     this->seed = seed;

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $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 +40 @@
 #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 +82 @@
 G4Abla::~G4Abla()
 {
+  delete fissionModel;
   delete randomGenerator;
   delete pace;
@@ -483 +457 @@
       //               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 +534 @@
       //               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 +1582 @@
                     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 +1732 @@
 
   imaxwell = 1;
-  inttype = 0;
   
   // limiting of excitation energy where fission occurs                    
@@ -2842 +2815 @@
   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 +3677 @@
 // 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();

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $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 +36 @@
 G4AblaFission::G4AblaFission()
 {
+  hazard = 0;
+  randomGenerator = 0;
 }
 
@@ -1086 +1088 @@
 }
 
-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();

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $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 +36 @@
 G4AblaFissionSimfis18::G4AblaFissionSimfis18()
 {
+  hazard = 0;
+  randomGenerator = 0;
 }
 
@@ -1412 +1414 @@
 }
 
-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();

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $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 +57 @@
   randomGenerator = new G4InclGeant4Random();
   //  randomGenerator = new G4Ranecu();
+  clearState();
 }
 
@@ -85 +86 @@
   //  randomGenerator = new G4Ranecu();
   randomGenerator = new G4InclGeant4Random();
+  clearState();
 }
 
@@ -149 +151 @@
 
   theLogger = 0;
+  clearState();
 }
 
 G4Incl::~G4Incl()
 {
+  delete be;
+  delete ps;
+  delete fermi;
+  delete qvp;
+
   delete randomGenerator;
   delete light_gaus_nuc;
@@ -385 +393 @@
 void G4Incl::processEventIncl(G4InclInput *input)
 {
+  //G4cout <<"Starting event " << eventnumber << G4endl;
   if(input == 0) {
     G4cerr <<"G4Incl fatal error: NULL pointer passed as input!" << G4endl;
@@ -393 +402 @@
   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 +525 @@
   }
 
-  if(calincl->bulletType() == 4) { // pi0
+  // pi0
+  if(calincl->bulletType() == 4) {
     mprojo = 134.9764;
     ap = 0.0;
@@ -414 +532 @@
   }
 
-  if(calincl->bulletType() == 5) { // pi-
+  // pi-
+  if(calincl->bulletType() == 5) {
     mprojo = 139.56995;
     ap = 0.0;
@@ -420 +539 @@
   }
 
-  // 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 +576 @@
   }
 
-  if(calincl->bulletType() == 9) { // Alpha
+  // Alpha
+  if(calincl->bulletType() == 9) {
     mprojo = 3727.42;
     ap = 4.0;
@@ -457 +583 @@
   }
 
-  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 +605 @@
   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 +623 @@
   }
 
-  /**
-   * 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 +680 @@
   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 +692 @@
 
   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 +740 @@
         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 +747 @@
         }
       }
-
+      
       if(kind[j] == 3) { 
         varntp->avv[j] = -1;
@@ -656 +827 @@
         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 +840 @@
 
       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 +885 @@
     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 +8331 @@
 }
 
-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 +8730 @@
     } // 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 +9024 @@
 }
 
+// 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 +9116 @@
 //       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]];

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $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 +116 @@
   incl->setInput(calincl);
 
-  G4InclInput::printProjectileTargetInfo(aTrack, theNucleus);
-  calincl->printInfo();
+  //  G4InclInput::printProjectileTargetInfo(aTrack, theNucleus);
+  //  calincl->printInfo();
 
 #ifdef DEBUGINCL
@@ -439 +439 @@
         }
       }
+      delete theFermiBreakupResult;
+      theFermiBreakupResult = 0;
+
       if(std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 
         G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl;

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $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 +116 @@
   vgsldin.close();
 
-  int A = 0, Zbegin = 0, Zend = 0;
   G4String str1, str2, str3;
   for(int i = 0; i < 500; i++) {
@@ -124 +123 @@
   }
   
+  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);
+      }
     }
   }

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $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 +103 @@
   }
 
-  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 +134 @@
       G4cout <<"Badly defined target after swapping. Falling back to normal (non-swapped) mode." << G4endl;
       calincl = new G4InclInput(aTrack, theNucleus, false);
+    }
     }
   } else {
@@ -489 +487 @@
         }
       }
+      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 +589 @@
         }
       }
+      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 +701 @@
   }
 
+  delete fermiBreakUp;
   delete calincl;
   calincl = 0;

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $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 +34 @@
 
 #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 +47 @@
 
   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 +64 @@
 G4InclCascadeInterface::~G4InclCascadeInterface()
 {
+  delete thePrecoModel;
+  delete theExcitationHandler;
+
   delete hazard;
   delete varntp;
-  delete inclInput;
   delete ws;
   delete mat;
@@ -68 +79 @@
 
   G4int particleI;
+  G4int bulletType = 0;
 
   // Print diagnostic messages: 0 = silent, 1 and 2 = verbose
@@ -91 +103 @@
   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 +117 @@
   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 +136 @@
   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 +150 @@
 
     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 +160 @@
         G4cout <<"G4InclCascadeInterface: Try number = " << nTries << G4endl; 
       }
-      incl->processEventIncl(inclInput);
+      incl->processEventIncl(calincl);
 
       if(verboseLevel > 1) {
@@ -151 +171 @@
        * 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 +196 @@
 
       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 +213 @@
 
 #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 +287 @@
       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 +325 @@
         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 +343 @@
         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 +363 @@
       }
     }
+
+      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 +484 @@
     }
 #endif
-    
+   
     varntp->ntrack = 0; // Clean up the number of generated particles in the event.
   }
@@ -379 +494 @@
     theResult.SetStatusChange(stopAndKill);
 
-    G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
     cascadeParticle = new G4DynamicParticle(theTableOfParticles->FindParticle(aTrack.GetDefinition()), aTrack.Get4Momentum());
 
@@ -388 +502 @@
     }
     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 +516 @@
     }
 
-    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 +544 @@
   }
 
+  delete calincl;
+  calincl = 0;
   return &theResult;
 } 

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $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 +31 @@
 // Aatos Heikkinen, HIP (project coordination)
 
+#include <vector>
+
 #include "G4InclLightIonInterface.hh"
+#include "G4FermiBreakUp.hh"
 #include "math.h"
 #include "G4GenericIon.hh"
@@ -39 +42 @@
 {
   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 +54 @@
   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 +82 @@
 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 +109 @@
   }
 
-  // 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 +161 @@
   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 +228 @@
     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 +249 @@
        * 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 +280 @@
       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 +312 @@
     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 +337 @@
           new G4DynamicParticle(G4Proton::ProtonDefinition(), momDirection, eKin);
         particleIdentified++;
+        baryonNumberBalanceInINCL -= 1;
+        chargeNumberBalanceInINCL -= 1;
       }
 
@@ -228 +345 @@
           new G4DynamicParticle(G4Neutron::NeutronDefinition(), momDirection, eKin);
         particleIdentified++;
+        baryonNumberBalanceInINCL -= 1;
       }
 
@@ -234 +352 @@
           new G4DynamicParticle(G4PionPlus::PionPlusDefinition(), momDirection, eKin);
         particleIdentified++;
+        chargeNumberBalanceInINCL -= 1;
       }
 
@@ -240 +359 @@
           new G4DynamicParticle(G4PionZero::PionZeroDefinition(), momDirection, eKin);
         particleIdentified++;
+        chargeNumberBalanceInINCL -= 0;
       }
 
@@ -246 +366 @@
           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 +393 @@
             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 +414 @@
     }
 
+    // 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 +622 @@
     cascadeParticle = new G4DynamicParticle(theTableOfParticles->FindParticle(aTrack.GetDefinition()), aTrack.Get4Momentum());
 
-    theResult.AddSecondary(cascadeParticle);
+    result.push_back(cascadeParticle);
 
     if(verboseLevel > 1) {
@@ -317 +635 @@
         G4cout <<"G4InclLightIonInterface: Unknown bullet type" << G4endl;      
         G4cout <<"Bullet particle name: " << cascadeParticle->GetDefinition()->GetParticleName() << G4endl;
-      }      
+      }
       if(verboseLevel > 3) {
         diagdata <<"G4InclLightIonInterface: Unknown bullet type" << G4endl;      
@@ -324 +642 @@
     }
 
-    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 +670 @@
   }
 
+  // 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;
 }