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

Timestamp:

Nov 5, 2010, 3:45:55 PM (15 years ago)

Author:

garnier

Message:

update ti head

Location:

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

Files:

• G4Abla.cc (modified) (6 diffs)
• G4AblaEvaporation.cc (modified) (4 diffs)
• G4AblaFissionBase.cc (added)
• G4Incl.cc (modified) (175 diffs)
• G4InclAblaCascadeInterface.cc (modified) (18 diffs)
• G4InclAblaLightIonInterface.cc (modified) (21 diffs)
• G4InclCascadeInterface.cc (modified) (13 diffs)
• G4InclLightIonInterface.cc (modified) (8 diffs)

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4Abla.cc,v 1.20 2009/11/18 10:43:14 kaitanie Exp $ 
+// $Id: G4Abla.cc,v 1.22 2010/10/26 02:47:59 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -98 +98 @@
 }
 
+void G4Abla::setVerboseLevel(G4int level)
+{
+  verboseLevel = level;
+  if(verboseLevel > G4InclUtils::silent) {
+    G4cout <<";; G4Abla: Setting verbose level to " << verboseLevel << G4endl;
+  }
+  fissionModel->setVerboseLevel(verboseLevel);
+}
+
 G4Abla::~G4Abla()
 {
@@ -111 +120 @@
 }
 
+void G4Abla::registerLogger(G4VInclLogger *theLogger) {
+  if(theLogger != NULL) {
+    this->theLogger = theLogger;
+  }
+}
+
 // Main interface to the evaporation
 
@@ -118 +133 @@
 // G4Fragment?
 
-void G4Abla::breakItUp(G4double nucleusA, G4double nucleusZ, G4double nucleusMass, G4double excitationEnergy,
+void G4Abla::breakItUp(G4int nucleusA, G4int nucleusZ, G4double nucleusMass, G4double excitationEnergy,
                        G4double angularMomentum, G4double recoilEnergy, G4double momX, G4double momY, G4double momZ,
                        G4int eventnumber)
@@ -187 +202 @@
   G4double esrem = excitationEnergy;
   
-  G4double aprf = nucleusA;
-  G4double zprf = nucleusZ;
+  G4double aprf = (double) nucleusA;
+  G4double zprf = (double) nucleusZ;
   G4double mcorem = nucleusMass;
   G4double ee = excitationEnergy;
@@ -208 +223 @@
   
   G4double pcorem = std::sqrt(erecrem*(erecrem +2.*938.2796*nucleusA));
+#ifdef G4INCLDEBUG
+  theLogger->fillHistogram1D("pcorem", pcorem);
+#endif
     // G4double pcorem = std::sqrt(std::pow(momX,2) + std::pow(momY,2) + std::pow(momZ,2));
   if(pcorem != 0) { // Guard against division by zero.

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4AblaEvaporation.cc,v 1.4 2008/10/24 21:07:40 dennis Exp $
+// $Id: G4AblaEvaporation.cc,v 1.6 2010/10/26 02:47:59 kaitanie Exp $
 //
 #include <numeric>
@@ -82 +82 @@
   hazard->igraine[17] = 33759;
   hazard->igraine[18] = 13227;
+
+  G4VarNtp *evaporationResult = new G4VarNtp();
+  G4Volant *volant = new G4Volant();
+
+  // Initialize evaporation.
+  abla = new G4Abla(hazard, volant, evaporationResult);
+  abla->initEvapora();
 }
 
@@ -108 +115 @@
 }
 
-G4FragmentVector * G4AblaEvaporation::BreakItUp(const G4Fragment &theNucleus) {
- 
-
-  G4VarNtp *varntp = new G4VarNtp();
-  G4Volant *volant = new G4Volant();
-
-  G4Abla *abla = new G4Abla(hazard, volant, varntp);
-  G4cout <<"Initializing evaporation..." << G4endl;
-  abla->initEvapora();
-  G4cout <<"Initialization complete!" << G4endl;
-  
-  G4double nucleusA = theNucleus.GetA();
-  G4double nucleusZ = theNucleus.GetZ();
+G4FragmentVector * G4AblaEvaporation::BreakItUp(const G4Fragment &theNucleus) {  
+  G4int nucleusA = theNucleus.GetA_asInt();
+  G4int nucleusZ = theNucleus.GetZ_asInt();
   G4double nucleusMass = G4NucleiProperties::GetNuclearMass(nucleusA, nucleusZ);
   G4double excitationEnergy = theNucleus.GetExcitationEnergy();
@@ -137 +134 @@
 
   varntp->ntrack = -1;
-  varntp->massini = theNucleus.GetA();
-  varntp->mzini = theNucleus.GetZ();
+  varntp->massini = theNucleus.GetA_asInt();
+  varntp->mzini = theNucleus.GetZ_asInt();
 
   std::vector<G4DynamicParticle*> cascadeParticles;

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4Incl.cc,v 1.30 2010/04/27 16:02:37 kaitanie Exp $ 
+// $Id: G4Incl.cc,v 1.33 2010/10/29 06:48:43 gunter Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -49 +49 @@
   densFunction = 5;
 
+  // Default: no Fermi break-up
+  useFermiBreakup = false;
+
+  // Default: no projectile spectators
+  useProjSpect = false;
+
   randomGenerator = new G4InclGeant4Random();
-  //randomGenerator = new G4Ranecu();
+  //  randomGenerator = new G4Ranecu();
 }
 
@@ -56 +62 @@
 {
   verboseLevel = 0;
+
+  // Default: no Fermi break-up
+  useFermiBreakup = false;
+
+  // Default: no projectile spectators
+  useProjSpect = false;
 
   // Set functions to be used for integration routine.  
@@ -71 +83 @@
   ws = aWs;
 
-  //randomGenerator = new G4Ranecu();
+  //  randomGenerator = new G4Ranecu();
   randomGenerator = new G4InclGeant4Random();
 }
 
-G4Incl::G4Incl(G4Hazard *aHazard, G4Calincl *aCalincl, G4Ws *aWs, G4Mat *aMat, G4VarNtp *aVarntp)
+G4Incl::G4Incl(G4Hazard *aHazard, G4InclInput *aCalincl, G4Ws *aWs, G4Mat *aMat, G4VarNtp *aVarntp)
 {
   verboseLevel = 0;
+
+  // Default: no Fermi break-up
+  useFermiBreakup = false;
+
+  // Default: no projectile spectators
+  useProjSpect = false;
 
   // Set functions to be used for integration routine.    
@@ -95 +113 @@
 
   randomGenerator = new G4InclGeant4Random();
-  // randomGenerator = new G4Ranecu();
+  //  randomGenerator = new G4Ranecu();
   light_gaus_nuc = new G4LightGausNuc();
   light_nuc = new G4LightNuc();
@@ -111 +129 @@
   bl10 = new G4Bl10();
   kindstruct = new G4Kind();
+  bev = new G4Bev();
   paul = new G4Paul();
   varavat = new G4VarAvat();
   varavat->kveux = 0;
+
+  be = new G4VBe();
+  ps = new G4InclProjSpect();
+  fermi = new G4InclFermi();
+  qvp = new G4QuadvectProjo();
 
   volant = new G4Volant();
@@ -123 +147 @@
   abla = new G4Abla(hazard, volant, evaporationResult);
   abla->initEvapora();
+
+  theLogger = 0;
 }
 
@@ -143 +169 @@
   delete bl10;
   delete kindstruct;
+  delete bev;
   delete paul;
   delete varavat;
@@ -149 +176 @@
   delete evaporationResult;
   delete volant;
+}
+
+void G4Incl::setUseFermiBreakUp(G4bool useIt)
+{
+  useFermiBreakup = useIt;
+}
+
+void G4Incl::setUseProjectileSpectators(G4bool useIt)
+{
+  useProjSpect = useIt;
 }
 
@@ -236 +273 @@
 {
   verboseLevel = level;
+  if(verboseLevel > G4InclUtils::silent) {
+    G4cout <<";; G4Incl: Setting verbose level to " << verboseLevel << G4endl;
+  }
+  abla->setVerboseLevel(level);
 }
 
@@ -268 +309 @@
 }
 
-void G4Incl::setCalinclData(G4Calincl *newCalincl)
+void G4Incl::setInput(G4InclInput *newCalincl)
 {
   calincl = newCalincl;
@@ -342 +383 @@
  */
 
-void G4Incl::processEventIncl()
-{
+void G4Incl::processEventIncl(G4InclInput *input)
+{
+  if(input == 0) {
+    G4cerr <<"G4Incl fatal error: NULL pointer passed as input!" << G4endl;
+    return;
+  }
+  calincl = input;
+
   const G4double uma = 931.4942;
   const G4double melec = 0.511;
@@ -355 +402 @@
 
   varntp->clear();
-  if(calincl->f[6] == 3.0) { // pi+ 
+  if(calincl->bulletType() == 3) { // pi+ 
     mprojo = 139.56995;
     ap = 0.0;
@@ -361 +408 @@
   }
 
-  if(calincl->f[6] == 4.0) { // pi0
+  if(calincl->bulletType() == 4) { // pi0
     mprojo = 134.9764;
     ap = 0.0;
@@ -367 +414 @@
   }
 
-  if(calincl->f[6] == 5.0) { // pi-
+  if(calincl->bulletType() == 5) { // pi-
     mprojo = 139.56995;
     ap = 0.0;
@@ -374 +421 @@
 
   // Coulomb en entree seulement pour les particules ci-dessous.
-  if(calincl->f[6] == 1.0) { // proton
+  if(calincl->bulletType() == 1) { // proton
     mprojo = 938.27231;
     ap = 1.0;
@@ -380 +427 @@
   }
   
-  if(calincl->f[6] == 2.0) { // neutron  
+  if(calincl->bulletType() == 2) { // neutron  
     mprojo = 939.56563;
     ap = 1.0;
@@ -386 +433 @@
   }
   
-  if(calincl->f[6] == 6.0) { // deuteron
+  if(calincl->bulletType() == 6) { // deuteron
     mprojo = 1875.61276;
     ap = 2.0;
@@ -392 +439 @@
   }
   
-  if(calincl->f[6] == 7.0) { // triton
+  if(calincl->bulletType() == 7) { // triton
     mprojo = 2808.95;
     ap = 3.0;
@@ -398 +445 @@
   }
   
-  if(calincl->f[6] == 8.0) { // He3
+  if(calincl->bulletType() == 8) { // He3
     mprojo = 2808.42;
     ap = 3.0;
@@ -404 +451 @@
   }
 
-  if(calincl->f[6] == 9.0) { // Alpha
+  if(calincl->bulletType() == 9) { // Alpha
     mprojo = 3727.42;
     ap = 4.0;
@@ -410 +457 @@
   }
 
-  pbeam = std::sqrt(calincl->f[2]*(calincl->f[2] + 2.0*mprojo));         
-
-  G4double at = calincl->f[0];
+  if(calincl->bulletType() == -12) { // Carbon
+    mprojo = 12.0 * 938.27231;
+    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();
+  }
+
+  pbeam = std::sqrt(calincl->bulletE()*(calincl->bulletE() + 2.0*mprojo));         
+
+  G4double at = calincl->targetA();
        
-  calincl->f[3] = 0.0;    // seuil sortie proton
-  calincl->f[7] = 0.0;    // seuil sortie neutron
-
   G4int ibert = 1;
 
@@ -435 +489 @@
   G4double probaTrans = 0.0;
   G4double rndm = 0.0;
-  if((calincl->f[6] == 1.0) || (calincl->f[6] >= 6.0)) {
-    probaTrans = coulombTransm(calincl->f[2],ap,zp,calincl->f[0],calincl->f[1]);
+
+  if((calincl->bulletType() == 1) || (calincl->bulletType() >= 6)) {
+    probaTrans = coulombTransm(calincl->bulletE(),ap,zp,calincl->targetA(),calincl->targetZ());
     standardRandom(&rndm, &(hazard->ial));
     if(rndm <= (1.0 - probaTrans)) {
@@ -447 +502 @@
    * Call the actual INCL routine.
    */
-  pnu(&ibert, &nopart,&izrem,&iarem,&esrem,&erecrem,&alrem,&berem,&garem,&bimpac,&jrem);
+
+  //  randomGenerator->printSeeds();
+  G4double jremx, jremy, jremz;
+  pnu(&ibert, &nopart,&izrem,&iarem,&esrem,&erecrem,&alrem,&berem,&garem,&bimpac,&jrem,
+      &jremx, &jremy, &jremz);
   forceAbsor(&nopart, &iarem, &izrem, &esrem, &erecrem, &alrem, &berem, &garem, &jrem);
   G4double aprf = double(iarem); // mass number of the prefragment
@@ -480 +539 @@
    * -following the approximations of the cugnon code (esrem...)
    */
-  G4double f0 = calincl->f[0];
-  G4double f1 = calincl->f[1];
-  G4double f2 = calincl->f[2];
+  G4double f0 = calincl->targetA();
+  G4double f1 = calincl->targetZ();
+  G4double f2 = calincl->bulletE();
   G4double mcorem = mprojo + f2 + abla->pace2(f0, f1) + f0 * uma - f1 * melec;
 
@@ -680 +739 @@
 
 
-void G4Incl::processEventInclAbla(G4int eventnumber)
-{
+void G4Incl::processEventInclAbla(G4InclInput *input, G4int eventnumber)
+{
+  //G4cout <<"Starting event " << eventnumber << G4endl;
+  if(input == 0) {
+    G4cerr <<"G4Incl fatal error: NULL pointer passed as input!" << G4endl;
+    return;
+  }
+  calincl = input;
+
   const G4double uma = 931.4942;
   const G4double melec = 0.511;
+  const G4double fmp = 938.2796;
 
   G4double pcorem = 0.0;
@@ -694 +761 @@
   varntp->clear();
 
+  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->f[6] == 3.0) { 
+  if(calincl->bulletType() == 3) { 
     mprojo = 139.56995;
     ap = 0.0;
@@ -702 +874 @@
 
   // pi0
-  if(calincl->f[6] == 4.0) {
+  if(calincl->bulletType() == 4) {
     mprojo = 134.9764;
     ap = 0.0;
@@ -709 +881 @@
 
   // pi-
-  if(calincl->f[6] == 5.0) {
+  if(calincl->bulletType() == 5) {
     mprojo = 139.56995;
     ap = 0.0;
@@ -718 +890 @@
 
   // proton
-  if(calincl->f[6] == 1.0) {
+  if(calincl->bulletType() == 1) {
     mprojo = 938.27231;
     ap = 1.0;
@@ -725 +897 @@
 
   // neutron  
-  if(calincl->f[6] == 2.0) {
+  if(calincl->bulletType() == 2) {
     mprojo = 939.56563;
     ap = 1.0;
@@ -732 +904 @@
 
   // deuteron
-  if(calincl->f[6] == 6.0) {
+  if(calincl->bulletType() == 6) {
     mprojo = 1875.61276;
     ap = 2.0;
@@ -739 +911 @@
 
   // triton
-  if(calincl->f[6] == 7.0) {
+  if(calincl->bulletType() == 7) {
     mprojo = 2808.95;
     ap = 3.0;
@@ -746 +918 @@
 
   // He3
-  if(calincl->f[6] == 8.0) {
+  if(calincl->bulletType() == 8) {
     mprojo = 2808.42;
     ap = 3.0;
@@ -753 +925 @@
 
   // Alpha
-  if(calincl->f[6] == 9.0) {
+  if(calincl->bulletType() == 9) {
     mprojo = 3727.42;
     ap = 4.0;
@@ -759 +931 @@
   }
 
-  pbeam = std::sqrt(calincl->f[2]*(calincl->f[2] + 2.0*mprojo));         
-
-  G4double at = calincl->f[0];
+  // Carbon
+  if(calincl->bulletType() == -12) {
+    mprojo = 6.0*938.27231 + 6.0*939.56563;
+    ap = 12.0;
+    zp = 6.0;
+  }
+
+  pbeam = std::sqrt(calincl->bulletE()*(calincl->bulletE() + 2.0*mprojo));         
+
+  G4double at = calincl->targetA();
        
-  calincl->f[3] = 0.0;    //     !seuil sortie proton
-  calincl->f[7] = 0.0;  //       !seuil sortie neutron
-
   G4int ibert = 1;
 
@@ -777 +953 @@
   G4double bimpac = 0.0;
   G4int jrem = 0;
+  G4double xjrem = 0.0, yjrem = 0.0, zjrem = 0.0;
   G4double alrem = 0.0;
 
@@ -784 +961 @@
   G4double rndm = 0.0;
 
-  if((calincl->f[6] == 1.0) || (calincl->f[6] >= 6.0)) {
-    //    probaTrans = coulombTransm(calincl->f[2],apro,zpro,calincl->f[0],calincl->f[1]);
-    probaTrans = coulombTransm(calincl->f[2],ap,zp,calincl->f[0],calincl->f[1]);
+  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));
     if(rndm <= (1.0 - probaTrans)) {
@@ -794 +971 @@
   }
 
+  //  G4cout <<"Before PNU:" << G4endl;
+  //  randomGenerator->printSeeds();
   // Call the actual INCL routine:
-  pnu(&ibert, &nopart,&izrem,&iarem,&esrem,&erecrem,&alrem,&berem,&garem,&bimpac,&jrem);
+  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;
@@ -846 +1032 @@
   //                -from energy balance (input - all emitted energies)
   //                -following the approximations of the cugnon code (esrem...)
-  G4double mcorem = mprojo + calincl->f[2] + abla->pace2(double(calincl->f[0]),double(calincl->f[1]))
-    + calincl->f[0]*uma - calincl->f[1]*melec;
+  G4double mcorem = mprojo + calincl->bulletE() + abla->pace2(double(calincl->targetA()),double(calincl->targetZ()))
+    + calincl->targetA()*uma - calincl->targetZ()*melec;
 
   G4double pxbil = 0.0;
@@ -854 +1040 @@
 
   if(nopart > -1) {
-    for(G4int j = 0; j < nopart; j++) {
-      if(ep[j] < 0.0) continue; // Workaround to avoid negative energies (and taking std::sqrt of a negative number).
-      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) { 
@@ -998 +1212 @@
     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) {
@@ -1013 +1234 @@
 
     // on recopie le remnant dans le ntuple
-    // varntp->ntrack = varntp->ntrack + 1;
+    varntp->ntrack = varntp->ntrack + 1;
     varntp->massini = iarem;
     varntp->mzini = izrem;
     varntp->exini = esrem;
+    varntp->pxrem = pxrem;
+    varntp->pyrem = pyrem;
+    varntp->pzrem = pzrem;
+    varntp->mcorem = mcorem;
+    varntp->erecrem = pcorem;
+    varntp->erecrem = erecrem;
 
     if(verboseLevel > 2) {
@@ -1022 +1249 @@
       varntp->dump();
     }
-    // Evaporation/fission:
-    evaporationResult->ntrack = 0;
-    abla->breakItUp(varntp->massini, varntp->mzini, mcorem, varntp->exini, varntp->jremn,
-                    erecrem, pxrem, pyrem, pzrem, eventnumber);
-
-    if(verboseLevel > 2) {
-      G4cout << __FILE__ << ":" << __LINE__ << "Dump evaporationResult after evaporation: " << G4endl;
-      evaporationResult->dump();
-    }
-    for(G4int evaporatedParticle = 0; evaporatedParticle < evaporationResult->ntrack; evaporatedParticle++) {
-      if(evaporationResult->avv[evaporatedParticle] == 0 && evaporationResult->zvv[evaporatedParticle] == 0) { //Fix: Skip "empty" particles with A = 0 and Z = 0
-        // continue;
-      }
-      varntp->kfis = evaporationResult->kfis;
-      varntp->itypcasc[varntp->ntrack] = evaporationResult->itypcasc[evaporatedParticle];
-      varntp->avv[varntp->ntrack] = evaporationResult->avv[evaporatedParticle];
-      varntp->zvv[varntp->ntrack]= evaporationResult->zvv[evaporatedParticle];
-      varntp->plab[varntp->ntrack] = evaporationResult->plab[evaporatedParticle];
-      varntp->enerj[varntp->ntrack] = evaporationResult->enerj[evaporatedParticle];
-      varntp->tetlab[varntp->ntrack] = evaporationResult->tetlab[evaporatedParticle];
-      varntp->philab[varntp->ntrack] = evaporationResult->philab[evaporatedParticle];
-      varntp->ntrack++;
-      if(verboseLevel > 3) {
-        G4cout <<"G4Incl: Returning evaporation result"           << G4endl;
-        G4cout <<"G4Incl: 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;
-      }
-    }
-    if(verboseLevel > 2) {
-      G4cout <<"G4Incl: ntrack = "          << varntp->ntrack << G4endl;
-      G4cout <<"G4Incl: Done extracting..." << G4endl;
-    }
+    // Maximum remnant for Geant4 Fermi break-up: A = 17, Z = 8
+    if(varntp->massini < 17 && varntp->exini < 1000 && varntp->massini > 0 && varntp->mzini > 0 && useFermiBreakup) { // Choose between Fermi break-up and ABLA
+      varntp->needsFermiBreakup = true;
+    } else {
+      varntp->needsFermiBreakup = false;
+      // Evaporation/fission:
+      evaporationResult->ntrack = 0;
+      abla->breakItUp(G4int(varntp->massini), G4int(varntp->mzini), mcorem, varntp->exini, varntp->jremn,
+                      erecrem, pxrem, pyrem, pzrem, eventnumber);
+
+      if(verboseLevel > 2) {
+        G4cout << __FILE__ << ":" << __LINE__ << "Dump evaporationResult after evaporation: " << G4endl;
+        evaporationResult->dump();
+      }
+      for(G4int evaporatedParticle = 0; evaporatedParticle < evaporationResult->ntrack; evaporatedParticle++) {
+        if(evaporationResult->avv[evaporatedParticle] == 0 && evaporationResult->zvv[evaporatedParticle] == 0) { //Fix: Skip "empty" particles with A = 0 and Z = 0
+          // continue;
+        }
+        varntp->kfis = evaporationResult->kfis;
+        varntp->itypcasc[varntp->ntrack] = evaporationResult->itypcasc[evaporatedParticle];
+        varntp->avv[varntp->ntrack] = evaporationResult->avv[evaporatedParticle];
+        varntp->zvv[varntp->ntrack]= evaporationResult->zvv[evaporatedParticle];
+        varntp->plab[varntp->ntrack] = evaporationResult->plab[evaporatedParticle];
+        varntp->enerj[varntp->ntrack] = evaporationResult->enerj[evaporatedParticle];
+        varntp->tetlab[varntp->ntrack] = evaporationResult->tetlab[evaporatedParticle];
+        varntp->philab[varntp->ntrack] = evaporationResult->philab[evaporatedParticle];
+        varntp->ntrack++;
+        if(verboseLevel > 3) {
+          G4cout <<"G4Incl: Returning evaporation result"           << G4endl;
+          G4cout <<"G4Incl: 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;
+        }
+      }
+      if(verboseLevel > 2) {
+        G4cout <<"G4Incl: ntrack = "          << varntp->ntrack << G4endl;
+        G4cout <<"G4Incl: Done extracting..." << G4endl;
+      }
+    }
+    if(calincl->bulletType() < 0 && useProjSpect && !useFermiBreakup) { // If we use projectile spectators for carbon but no fermi breakup
+      // Evaporation/fission:
+      evaporationResult->ntrack = 0;
+      //      G4cout <<"Warning: Using ABLA to de-excite projectile spectator..." << G4endl;
+      abla->breakItUp(G4int(varntp->masp), G4int(varntp->mzsp), ps->m_projspec, varntp->exsp, 0.0,
+                      ps->t_projspec, ps->p1_projspec, ps->p1_projspec, ps->p1_projspec, eventnumber);
+
+      if(verboseLevel > 2) {
+        G4cout << __FILE__ << ":" << __LINE__ << "Dump evaporationResult after evaporation: " << G4endl;
+        evaporationResult->dump();
+      }
+      for(G4int evaporatedParticle = 0; evaporatedParticle < evaporationResult->ntrack; evaporatedParticle++) {
+        if(evaporationResult->avv[evaporatedParticle] == 0 && evaporationResult->zvv[evaporatedParticle] == 0) { //Fix: Skip "empty" particles with A = 0 and Z = 0
+        }
+        varntp->kfis = evaporationResult->kfis;
+        varntp->itypcasc[varntp->ntrack] = evaporationResult->itypcasc[evaporatedParticle];
+        varntp->avv[varntp->ntrack] = evaporationResult->avv[evaporatedParticle];
+        varntp->zvv[varntp->ntrack]= evaporationResult->zvv[evaporatedParticle];
+        varntp->plab[varntp->ntrack] = evaporationResult->plab[evaporatedParticle];
+        varntp->enerj[varntp->ntrack] = evaporationResult->enerj[evaporatedParticle];
+        varntp->tetlab[varntp->ntrack] = evaporationResult->tetlab[evaporatedParticle];
+        varntp->philab[varntp->ntrack] = evaporationResult->philab[evaporatedParticle];
+        varntp->ntrack++;
+        if(verboseLevel > 3) {
+          G4cout <<"G4Incl: Returning evaporation result"           << G4endl;
+          G4cout <<"G4Incl: 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;
+        }
+      }
+    }
+#ifdef G4INCLDEBUG
+    theLogger->fillHistogram1D("bimpact", varntp->bimpact);
+#endif
+
+#ifdef G4INCLDEBUG
+    theLogger->fillHistogram1D("mzini", varntp->mzini);
+#endif
   }
   if(nopart == -2) {
@@ -1403 +1676 @@
     
   if(tz < 0) {
+#ifdef G4INCLDEBUG
+  theLogger->fillHistogram1D("interpolationResult", (saxw->y[0][saxw->imat] + saxw->s[0][saxw->imat]*tz));
+  theLogger->fillHistogram2D("interpolationPoints", (saxw->y[0][saxw->imat] + saxw->s[0][saxw->imat]*tz), xv);
+#endif
     return (saxw->y[0][saxw->imat] + saxw->s[0][saxw->imat]*tz);
   }
   else if(tz == 0) {
-    return (saxw->y[0][saxw->imat]);
+#ifdef G4INCLDEBUG
+  theLogger->fillHistogram1D("interpolationResult", (saxw->y[0][saxw->imat]));
+  theLogger->fillHistogram2D("interpolationPoints", (saxw->y[0][saxw->imat]), xv);
+#endif
+  return (saxw->y[0][saxw->imat]);
   }
   else { // tz > 0
@@ -1417 +1698 @@
     }
     if(tz >= 0) {
+#ifdef G4INCLDEBUG
+      theLogger->fillHistogram1D("interpolationResult", saxw->y[j][saxw->imat]);
+      theLogger->fillHistogram2D("interpolationPoints", saxw->y[j][saxw->imat], xv);
+#endif
       return saxw->y[j][saxw->imat];
     } else if(tz < 0.0) {
       j = j - 1;
       G4double dgx = xv - saxw->x[j][saxw->imat];
+#ifdef G4INCLDEBUG
+      theLogger->fillHistogram1D("interpolationResult", (saxw->y[j][saxw->imat] + saxw->s[j][saxw->imat]*dgx));
+      theLogger->fillHistogram2D("interpolationPoints", (saxw->y[j][saxw->imat] + saxw->s[j][saxw->imat]*dgx), xv);
+#endif
       return(saxw->y[j][saxw->imat] + saxw->s[j][saxw->imat]*dgx);
     }
@@ -1673 +1962 @@
 void G4Incl::pnu(G4int *ibert_p, G4int *nopart_p, G4int *izrem_p, G4int *iarem_p, G4double *esrem_p,
                  G4double *erecrem_p, G4double *alrem_p, G4double *berem_p, G4double *garem_p,
-                 G4double *bimpact_p, G4int *l_p)
-{
+                 G4double *bimpact_p, G4int *l_p, G4double *xjrem, G4double *yjrem, G4double *zjrem)
+{
+  //  G4cout <<"Now in pnu..." << G4endl;
+  //  G4cout <<"(clean-up)" << G4endl;
+  G4int npenter = 0;
+  G4int nnenter = 0;
+  G4int n_enter_pot = 0;
+  G4int avatarCounter = 0;
+  ps->clear(); // For projectile spectators
+  //  G4int i_c = 0;
+  G4double p1spec=0.;
+  G4double p2spec=0.;
+  G4double p3spec=0.;
+  //  G4double p_spec2 = 0.0;
+  //  G4double s_spec = 0.0;
+  // G4double e_spec = 0.0;
+  G4double xl1_i, xl2_i, xl3_i;
+  G4int idq = 0;
+  G4double destar = 0.0;
+  //  G4cout <<"(checkpoint 'setipszero0)" << G4endl;
+  G4int ips = 0;
+  //  G4double p1_spec = 0.0;
+  //  G4double p2_spec = 0.0;
+  //  G4double p3_spec = 0.0;
+  G4int n_activnuc = 0;
   G4int ibert = (*ibert_p);
   //  float f[15]; // = (*f_p);
@@ -1702 +2014 @@
   G4double b3 = 0.0;
   G4double bb2 = 0.0;
-  G4double be = 0.0;
+  G4double bevar = 0.0;
   G4double bmass[2000]; 
   for(G4int init_i = 0; init_i < 2000; init_i++) {
@@ -1736 +2048 @@
   G4double ener_max = 0.0;
   G4double eout = 0.0;
-  G4double eps_c[BL1SIZE]; 
-  for(G4int init_i = 0; init_i < BL1SIZE; init_i++) {
-    eps_c[init_i] = 0.0;
-  }
+  //  G4double eps_c[BL1SIZE]; 
+  //  G4double p3_c[BL1SIZE];
+  //  for(G4int init_i = 0; init_i < BL1SIZE; init_i++) {
+    //    eps_c[init_i] = 0.0;
+    //    p3_c[init_i] = 0.0;
+  //  }
   G4double epsv = 0.0;
   G4double erecg = 0.0;
@@ -1778 +2092 @@
   G4int imin = 0;
   G4int indic[3000]; 
+  G4int nb_transprojo=0;
+  G4double v_proj = 0.0;
   for(G4int init_i = 0; init_i < 3000; init_i++) {
     indic[init_i] = 0;
@@ -1825 +2141 @@
   G4int nbtest = 0;
   G4int nc[300];
+  G4bool isPartOfSpectatorNucleus[300];
   G4int npproj[300];
   for(G4int init_i = 0; init_i < 300; init_i++) {
@@ -1942 +2259 @@
   G4double xleng = 0.0;
   G4double xlengm = 0.0;
-  G4double xmodp = 0.0;
   G4double xpb = 0.0;
   G4double xq = 0.0;
@@ -1960 +2276 @@
   G4double xye = 0.0;
   G4double y = 0.0;
-  G4double p3_c[BL1SIZE];
   G4double q1[BL1SIZE];
   G4double q2[BL1SIZE];
@@ -2220 +2535 @@
    for(G4int i = 0; i < 300; i++) {
      npproj[i] = 0;
+     isPartOfSpectatorNucleus[i] = false;
      nc[i] = 0;
    }
@@ -2270 +2586 @@
   //    427     c     k6=0 no angular momentum conservation                             p-n02070
   //    428     c     k6=1 angular momentum conservation                                p-n02080
+  //                  Kclst=F(10) Light clusters computed if not =0
   //    429     c                                                                       p-n02090
   //    430     c     b=impact parameter                                                p-n02100
@@ -2312 +2629 @@
   // end for logging
 
+  G4int ia_breakup = 0;
   G4int jparticip[BL1SIZE];
   for(G4int i = 0; i < BL1SIZE; i++) {
@@ -2318 +2636 @@
 
   G4double beproj = 0.;
-  bl3->ia2 = G4int(std::floor(calincl->f[0] + 0.1)); // f(1)->f[0] and so on..., calincl added
-  G4int iz2 = G4int(std::floor(calincl->f[1] + 0.1));
+  bl3->ia2 = calincl->targetA(); // f(1)->f[0] and so on..., calincl added
+  G4int iz2 = calincl->targetZ();
   G4double r02 = 1.12;
-  kindstruct->kindf7 = int(std::floor(calincl->f[6] + 0.1));
-
-  bl3->ia1 = ia1t[G4int(kindstruct->kindf7)-1];
-  G4int iz1 = iz1t[G4int(kindstruct->kindf7)-1];
-  G4double fmpinc = fmpinct[G4int(kindstruct->kindf7)-1];
-  G4double rms1 = light_gaus_nuc->rms1t[G4int(kindstruct->kindf7)-1];
-  G4double pf1 = light_gaus_nuc->pf1t[G4int(kindstruct->kindf7)-1];
-  G4double tlab = calincl->f[2];
+  kindstruct->kindf7 = calincl->bulletType();
+  if(kindstruct->kindf7 < 0) kindstruct->kindf7 = kindstruct->kindf7-1;
+  G4int kclst = calincl->getClusterOption();
+
+  // G4cout <<"Projectile kind = " << kindstruct->kindf7 << G4endl;
+
+  G4int iz1 = 0;
+  G4double fmpinc = 0.0;
+  G4double rms1 = 0.0;
+  G4double pf1 = 0.0;
+  G4double tlab = calincl->bulletE();;
+
+  if(kindstruct->kindf7 > 0) {
+    bl3->ia1 = ia1t[G4int(kindstruct->kindf7)-1]; // In the C++ version indices start from 0
+    iz1 = iz1t[G4int(kindstruct->kindf7)-1];
+    fmpinc = fmpinct[G4int(kindstruct->kindf7)-1];
+    rms1 = light_gaus_nuc->rms1t[G4int(kindstruct->kindf7)-1];
+    pf1 = light_gaus_nuc->pf1t[G4int(kindstruct->kindf7)-1];
+
+  } else {
+   // Starting values for carbon beams:
+    if(kindstruct->kindf7 == -12) { // Handle Carbon-12 as a special case
+      bl3->ia1 = std::abs(kindstruct->kindf7 + 1);
+      // check here that kindf7 is -a and ia1 is a of the projectile.
+      //                        write(6,*) 'kindf7,ia1:',kindf7,ia1
+      be->ia_be=bl3->ia1;
+      iz1=be->iz_be;
+      fmpinc=bl3->ia1*fmp - be->bind_be;
+      rms1=be->rms_be;
+      pf1=be->pms_be;
+      ia_breakup = 10; // Stopping threshold for the cascade in case of light nuclei
+      //      ia_breakup=calincl->getCascadeStoppingAThreshold();
+    } else { // For extended projectiles
+      bl3->ia1 = calincl->extendedProjectileA();
+      be->ia_be = calincl->extendedProjectileA();
+      iz1=be->iz_be;
+      fmpinc=bl3->ia1*fmp - be->bind_be;
+      rms1=be->rms_be;
+      pf1=be->pms_be;
+      ia_breakup = 10; // Stopping threshold for the cascade in case of light nuclei
+      //      ia_breakup=calincl->getCascadeStoppingAThreshold();
+    }
+  }
+
+  if(verboseLevel > 1) {
+    G4cout <<"rms1 = " << rms1 << G4endl;
+    G4cout <<"pf1 = " << pf1 << G4endl;
+  }
 
   G4int k2 = 0;
@@ -2338 +2696 @@
 
   // material number:      
-  saxw->imat = G4int(std::floor(calincl->f[8] + 0.5)); // f(9) -> f[8]
+  saxw->imat = 0;
   // espace de phases test (r et p) pour pauli: 
   // valeur recommandee par j.c. v-test=0.589 h**3:
@@ -2356 +2714 @@
   // temfin (time at which the inc is stopped), tmax5 defined after chosing b
 
-  G4double v0 = calincl->f[4]; // f(5)->f[4]
+  G4double v0 = calincl->getPotential();
   G4double v1 = v0;
   bl8->rathr = 0.;
@@ -2428 +2786 @@
     G4cout <<"Radius bl3->r2 = " << bl3->r2 << G4endl;
   }
-  
+
   G4double tnr = tlab;
   G4double binc = std::sqrt(tnr*tnr + 2.0*tlab*fmpinc)/(tnr+fmpinc);
-  G4double ginc=1.0/std::sqrt(1.0 - binc*binc);
+  G4double ginc = 1.0/std::sqrt(1.0 - binc*binc);
   G4double pinc = fmpinc*binc*ginc;
 
@@ -2463 +2821 @@
   //    710     c      if (kindf7.gt.2) bmax=bmax       ! a.b. (avec w.s., idem les nucleons)
   // deutons cv 22/01/2001
-  if (kindstruct->kindf7 <= 2)  { //then
+  if (kindstruct->kindf7 < 6 && kindstruct->kindf7 > 0)  { //then
     ws->bmax = ws->bmax;     // comme alain
   }
   else {
-    if (kindstruct->kindf7 < 6) { // then    
-      ws->bmax = ws->bmax;   // comme alain
-    }
-    else {
-      beproj = fmpinc - bl3->ia1*fmp;
-      ws->bmax = ws->rmaxws + rms1;     // maximum extension of the nucleus ( w.s.)
-    }
+    beproj = fmpinc - bl3->ia1*fmp;
+    ws->bmax = ws->rmaxws + rms1;     // maximum extension of the nucleus ( w.s.)
   }
 
@@ -2496 +2849 @@
       bred = 0.;
     }
-    if (kindstruct->kindf7 <= 2) {
+    if (kindstruct->kindf7 <= 2 && kindstruct->kindf7 > 0) {
       if (tlab < 400.0) {
         cb0 = 6.86 - 0.0035 * tlab;
@@ -2515 +2868 @@
     }
     else {
-      if (kindstruct->kindf7 < 6) {
+      if (kindstruct->kindf7 < 6 && kindstruct->kindf7 > 0) {
         // here for pions:
         temfin = 30.0*std::pow((float(bl3->ia2)/208.0),0.25)*(1.0 - 0.2*bred)*(1.0 - tlab/1250.0);
@@ -2557 +2910 @@
       // temfin=60.*(float(ia2)/208.)**0.25*(1.-0.2*bred)*(1.-tlab/1250.)
       // modified in april 2003 (more reasonable but not yet checked!)
-      temfin = 25.5*std::pow(G4double(bl3->ia2),0.16);  // pb208->60fm/c
+      //      temfin = 25.5*std::pow(G4double(bl3->ia2),0.16);  // pb208->60fm/c
+
+      // c modified in June 2005 (function of A and TLAB)
+      temfin = 30.18*std::pow(bl3->ia2,0.17)*(1.0-5.7*std::pow(10.0,(-5))*tlab);
     }
     else {
@@ -2567 +2923 @@
   // deutons
   // a way to change stopping time f[5] not used here
-  factemp = calincl->f[5]; // f(6)->f[5]
+  factemp = calincl->getTimeScale();
   // attention !!! 30/04/2001 scaling time is now a multiplication factor
   temfin = temfin*factemp;
@@ -2604 +2960 @@
   else {
     npion = 1;
-    ipi[1] = 8 - 2*(kindstruct->kindf7);
+    if(kindstruct->kindf7 == -666) { // Extended projectile...
+      ipi[1] = 8 - 2*be->ia_be;
+    } else {
+      ipi[1] = 8 - 2*(kindstruct->kindf7);
+    }
     y1[1] = 0.0;
     y2[1] = 0.0;
@@ -2617 +2977 @@
   goto pnu9;
   //    850     
- pnu7: 
-  s1t1 = 0.0;
-  s2t1 = 0.0;
-  s3t1 = 0.0;
-  sp1t1 = 0.0;
-  sp2t1 = 0.0;
-  sp3t1 = 0.0;
-  for(G4int i = 1; i <= bl3->ia1; i++) {
-  //for(G4int i = 1; i <= (bl3->ia1-1); i++) {
-  //  for(G4int i = 1; i < (bl3->ia1-1); i++) {
-    bl9->hel[i] = 0;
-    bl5->nesc[i] = 0;
-    bl1->ind2[i] = 1;
-    if (i > iz1) {
-      bl1->ind2[i] = -1;
-    }
-    bl1->ind1[i] = 0;
-    // deutons
-    jparticip[i] = 1;
-
-    // deutons 
+ pnu7:
+  //  G4cout <<"(checkpoint 'pnu7)" << G4endl;
+  if(kindstruct->kindf7 == 6) { // Deuteron
     gaussianRandom(&xga);
-    bl3->x1[i] = xga*rms1*0.57735;
-    s1t1 = s1t1 + bl3->x1[i];
+    bl3->x1[1] = xga*rms1*0.5775;
     gaussianRandom(&xga);
-    bl3->x2[i] = xga*rms1*0.57735;
-    s2t1 = s2t1 + bl3->x2[i];
+    bl3->x2[1] = xga*rms1*0.5775;
     gaussianRandom(&xga);
-    bl3->x3[i] = xga*rms1*0.57735; 
-    s3t1 = s3t1 + bl3->x3[i];
-
-    if(kindstruct->kindf7 == 6) { //then
-      // deuteron density from paris potential in q space:      
-      standardRandom(&xq, &(hazard->igraine[9]));
-      qdeut = splineab(xq) * 197.3289;
-      standardRandom(&u, &(hazard->igraine[10]));
-      cstet = u*2 - 1;
-      sitet = std::sqrt(1.0 - std::pow(cstet,2));
-      standardRandom(&v, &(hazard->igraine[11]));
-      phi = 2.0*3.141592654*v;
-
-      bl1->p1[i] = qdeut*sitet*std::cos(phi);
-      bl1->p2[i] = qdeut*sitet*std::sin(phi);
-      bl1->p3[i] = qdeut*cstet;
-    }
-    else {
-      // density of composite as a gaussien in q space:
+    bl3->x3[1] = xga*rms1*0.5775;
+
+    bl3->x1[2] = -bl3->x1[1];
+    bl3->x2[2] = -bl3->x2[1];
+    bl3->x3[2] = -bl3->x3[1];
+
+    // Deuteron density from Paris potential in q space:
+    standardRandom(&xq, &(hazard->igraine[9]));
+    qdeut = splineab(xq) * 197.3289;
+    standardRandom(&u, &(hazard->igraine[11]));
+    cstet = u * 2 - 1;
+    sitet = std::sqrt(1.0 - std::pow(cstet, 2));
+    standardRandom(&v, &(hazard->igraine[11]));
+    phi = 2.0 * 3.141592654 * v;
+    bl1->p1[1] = qdeut * sitet * std::cos(phi);
+    bl1->p2[1] = qdeut * sitet * std::sin(phi);
+    bl1->p3[1] = qdeut * cstet;
+    bl1->eps[1] = w(bl1->p1[1], bl1->p2[1], bl1->p3[1], fmp);
+    bl1->p1[2] = - bl1->p1[1];
+    bl1->p2[2] = - bl1->p2[1];
+    bl1->p3[2] = - bl1->p3[1]; 
+    bl1->eps[2] = bl1->eps[1];
+    //    bl1->eps[2] = w(bl1->p1[2], bl1->p2[2], bl1->p3[2], fmp);
+
+    jparticip[1] = 1;
+    bl9->hel[1] = 0;
+    bl5->nesc[1] = 0;
+    bl1->ind2[1] = 1;
+    bl1->ind1[1] = 0;
+    jparticip[2] = 1;
+    bl9->hel[2] = 0;
+    bl5->nesc[2] = 0;
+    bl1->ind2[2] = -1;
+    bl1->ind1[2] = 0;
+  } else {
+    // Composite heavier than a deuteron
+    s1t1 = 0.0;
+    s2t1 = 0.0;
+    s3t1 = 0.0;
+    sp1t1 = 0.0;
+    sp2t1 = 0.0;
+    sp3t1 = 0.0;
+    G4double sei = 0.0;
+
+    bl9->hel[0] = 0;
+    bl5->nesc[0] = 0;
+    for(G4int i = 1; i <= bl3->ia1; i++) {
+      bl9->hel[i] = 0;
+      bl5->nesc[i] = 0;
+      bl1->ind2[i] = 1;
+      if(i > iz1) {
+        bl1->ind2[i] = -1;
+      }
+      bl1->ind1[i] = 0;
+      jparticip[i] = 1;
       gaussianRandom(&xga);
-      bl1->p1[i] = xga*pf1*0.57735;
+      bl3->x1[i] = xga * rms1 * 0.57735;
+      s1t1 = s1t1 + bl3->x1[i];
       gaussianRandom(&xga);
-      bl1->p2[i] = xga*pf1*0.57735;
+      bl3->x2[i] = xga * rms1 * 0.57735;
+      s2t1 = s2t1 + bl3->x2[i];
       gaussianRandom(&xga);
-      bl1->p3[i] = xga*pf1*0.57735;
-    }
-
-    bl1->eps[i] = w(bl1->p1[i],bl1->p2[i],bl1->p3[i],fmp);
-
-    sp1t1 = sp1t1 + bl1->p1[i];
-    sp2t1 = sp2t1 + bl1->p2[i];
-    sp3t1 = sp3t1 + bl1->p3[i];
-  }
-
-  bl9->hel[bl3->ia1] = 0;
-  bl5->nesc[bl3->ia1] = 0;
-  bl1->ind2[bl3->ia1] = -1;
-  bl1->ind1[bl3->ia1] = 0;
-  bl3->x1[bl3->ia1] = -s1t1;
-  bl3->x2[bl3->ia1] = -s2t1;
-  bl3->x3[bl3->ia1] = -s3t1;
-  bl1->p1[bl3->ia1] = -sp1t1;
-  bl1->p2[bl3->ia1] = -sp2t1;
-  bl1->p3[bl3->ia1] = -sp3t1;
-  bl1->eps[bl3->ia1] = w(bl1->p1[bl3->ia1],bl1->p2[bl3->ia1],bl1->p3[bl3->ia1],fmp);
-  
-  // deutons
-  jparticip[bl3->ia1] = 1;
-  if(verboseLevel > 3) {
-    G4cout <<"Particle " << bl3->ia1-1 << " is now participant." << G4endl;
-  }
+      bl3->x3[i] = xga * rms1 * 0.57735;
+      s3t1 = s3t1 + bl3->x3[i];
+
+      G4int igaussfermi = 2;
+      if(igaussfermi == 1) { // Gaussian
+        // Density of composite as a gaussian in q
+        // space:
+        gaussianRandom(&xga);
+        bl1->p1[i] = xga * pf1 * 0.57735;
+        gaussianRandom(&xga);
+        bl1->p2[i] = xga * pf1 * 0.57735;
+        gaussianRandom(&xga);
+        bl1->p3[i] = xga * pf1 * 0.57735;
+      } else { // Fermi
+        // Density of composite as a Fermi sphere in q space:
+        for(int iloc=1; iloc <= 3; iloc++) {
+          standardRandom(&t[iloc],&(hazard->igraine[11]));
+        }
+        t[2]=-1+2.*t[2];
+        t[3]=6.283185*t[3];
+        t2=t[2];
+        t3=std::sqrt(1.-t2*t2);
+        t4=std::cos(t[3]);
+        t5=std::sin(t[3]);
+        y=fermi->pf*std::pow(t[1],0.33333333);
+        bl1->p1[i]=y*t3*t4;
+        bl1->p2[i]=y*t3*t5;
+        bl1->p3[i]=y*t2;
+      }
+      bl1->eps[i] = w(bl1->p1[i], bl1->p2[i], bl1->p3[i], fmp);
+      sei = sei + bl1->eps[i];
+
+      if(verboseLevel > 1) {
+        G4cout <<"projectile nucleon = " << i << G4endl;
+        G4cout <<"x1 = " << bl3->x1[i] << G4endl;
+        G4cout <<"x2 = " << bl3->x2[i] << G4endl;
+        G4cout <<"x3 = " << bl3->x3[i] << G4endl;
+        G4cout <<"p1 = " << bl1->p1[i] << G4endl;
+        G4cout <<"p2 = " << bl1->p2[i] << G4endl;
+        G4cout <<"p3 = " << bl1->p3[i] << G4endl;
+        G4cout <<"eps = " << bl1->eps[i] << G4endl;
+      }
+      sp1t1 = sp1t1 + bl1->p1[i];
+      sp2t1 = sp2t1 + bl1->p2[i];
+      sp3t1 = sp3t1 + bl1->p3[i];
+    }
+
+    s1t1 = s1t1/double(bl3->ia1);
+    s2t1 = s2t1/double(bl3->ia1);
+    s3t1 = s3t1/double(bl3->ia1);
+    sp1t1 = sp1t1/double(bl3->ia1);
+    sp2t1 = sp2t1/double(bl3->ia1);
+    sp3t1 = sp3t1/double(bl3->ia1);
+    sei = 0.0;
+
+    for(G4int i = 1; i <= bl3->ia1; i++) {
+      bl3->x1[i] = bl3->x1[i] - s1t1;
+      bl3->x2[i] = bl3->x2[i] - s2t1;
+      bl3->x3[i] = bl3->x3[i] - s3t1;
+      bl1->p1[i] = bl1->p1[i] - sp1t1;
+      bl1->p2[i] = bl1->p2[i] - sp2t1;
+      bl1->p3[i] = bl1->p3[i] - sp3t1;
+      bl1->eps[i] = w(bl1->p1[i], bl1->p2[i], bl1->p3[i], fmp);
+      sei = sei + bl1->eps[i];
+    }
+    v_proj = 0.0;
+    if(bl3->ia1 > 4) v_proj = (sei - fmpinc)/bl3->ia1;
+
+    //    bl1->dump(26);
+    //    bl3->dump();
+    //    G4cout <<"(checkpoint 'endheavyprojinit)" << G4endl;
+  } // End of: Composites heavier than deuteron
+
+
  pnu9: // continue
+  //  G4cout <<"(checkpoint 'pnu9)" << G4endl;
   // deutons
   // target preparation for 1 < a < 5 (with sum of momentum =0)
@@ -2860 +3286 @@
   // random azimuthal direction of the impact parameter (sept 99)
 
-  if (kindstruct->kindf7 <= 2) {
-    G4long testseed = hazard->igraine[13];
-    //    standardRandom(&tbid, &(hazard->igraine[13]));
-    standardRandom(&tbid, &testseed);
-    hazard->igraine[13] = testseed;
+  if (kindstruct->kindf7 <= 2 && kindstruct->kindf7 > 0) {
+    standardRandom(&tbid, &(hazard->igraine[13]));
     tbid = tbid*6.283185;
     bl3->x1[1] = bl3->x1[1] + b*std::cos(tbid);  //x1(1)->x1[1]                                        
     bl3->x2[1] = bl3->x2[1] + b*std::sin(tbid); //x2(1)->x2[1]
     bl3->x3[1] = bl3->x3[1] - z;
+    // Counter of p(n) entering the potential 4/2008 AB
+    if(bl1->ind2[1] == 1) { // then
+      npenter=1;
+      nnenter=0;
+    } else {
+      npenter=0;
+      nnenter=1;
+    } // endif
+
     // pour le ntuple des avatars:
     if(varavat->kveux == 1) {
@@ -2877 +3309 @@
   }
   else {
-    if (kindstruct->kindf7 < 6) { //then ! pour les pions on laisse
+    if (kindstruct->kindf7 < 6 && kindstruct->kindf7 > 0) { //then ! pour les pions on laisse
       //call standardRandom(tbid,iy14)
       standardRandom(&tbid, &(hazard->igraine[13]));
@@ -2916 +3348 @@
       //    pnu21:
       if (nmiss == bl3->ia1) { //then
+        //      G4cout <<"(checkpoint 'projectilemissedtarget)" << G4endl;
         nopart = -1;
         //      return;
@@ -2925 +3358 @@
       }
       else {
+        // Counter of p(n) entering the potential 4/2008 AB
+        if(bl1->ind2[ilm]==1) { //then
+          npenter=1;
+          nnenter=0;
+        } else {
+          npenter=0;
+          nnenter=1;
+        } // endif
+
         zshif = bl3->x3[ilm] - ztouch;
         standardRandom(&tbid, &(hazard->igraine[13]));
@@ -2943 +3385 @@
   }
 
+  // for rho(r),rho(q) checking
+#ifdef G4INCLDEBUG
+  for(G4int rho_i = bl3->ia1 + 1; rho_i <= ia; rho_i++) { // DO I=IA1+1,IA
+    G4double r_dist = std::sqrt(bl3->x1[rho_i] * bl3->x1[rho_i] + bl3->x2[rho_i] * bl3->x2[rho_i] + bl3->x3[rho_i] * bl3->x3[rho_i]);
+    G4double p_mod = std::sqrt(std::pow(bl1->p1[rho_i],2)+std::pow(bl1->p2[rho_i],2) + std::pow(bl1->p3[rho_i],2));
+    theLogger->fillHistogram1D("r_distrib", r_dist);
+    theLogger->fillHistogram1D("p_distrib", p_mod);
+    theLogger->fillHistogram2D("r-p_correl", r_dist, p_mod);
+  }
+#endif
+
   // initial momentum for all type of incident particles:
   xl1 = b*pinc*std::sin(tbid);                                  
@@ -2951 +3404 @@
   // (here,=lab frame)
 
-  be = 0.0;
+  bevar = 0.0;
   ge = 1.0;
-  b1 = (binc - be)/(1.0 - be*binc);
-  b2 = -be;
+  b1 = (binc - bevar)/(1.0 - bevar*binc);
+  b2 = -bevar;
   g1 = 1.0/std::sqrt(1.0 - b1*b1);
   g2 = 1.0;
   // deutons
   // here for nucleons
-  if (kindstruct->kindf7 <= 2) {
+  if (kindstruct->kindf7 <= 2 && kindstruct->kindf7 > 0) {
     bl1->eps[1] = g1*fmp + v0;
     bl1->p3[1] = std::sqrt(std::pow(bl1->eps[1],2) - std::pow(fmp,2));
@@ -2965 +3418 @@
   else {
     // here for pions
-    if (kindstruct->kindf7 < 6) { //then
+    if (kindstruct->kindf7 < 6 && kindstruct->kindf7 > 0) { //then
       q4[1] = g1*fmpi; // q4(1)->q4[0]
       q3[1] = b1*q4[1];
     }
-    else {
+    else { // Composite projectiles
       // here for composite projectiles:
       // the kinetic energy is below the threshold. put all
@@ -2985 +3438 @@
       // here the composit is above threshold
       for(G4int i = 1; i <= bl3->ia1; i++) { //do i=1,bl3->ia1  !save e,p in the composit rest frame
-        eps_c[i] = bl1->eps[i];
-        p3_c[i] = bl1->p3[i];
+        qvp->t_c[i] = bl1->eps[i] - fmp;
+        bl1->eps[i] = bl1->eps[i] - v_proj;
+        qvp->eps_c[i] = bl1->eps[i];
+        qvp->p3_c[i] = bl1->p3[i];
       } //enddo
 
@@ -2995 +3450 @@
       iflag = 0;
     pnu1870:
+      //      G4cout <<"(checkpoint 'pnu1870)" << G4endl;
       sueps = 0.0;
 
@@ -3000 +3456 @@
 
       for(G4int i = 1; i <= bl3->ia1; i++) { //do i=1,bl3->ia1
-        tte = eps_c[i];
-        bl1->eps[i] = g1*(eps_c[i] + b1*p3_c[i]);
-        bl1->p3[i] = g1*(b1*tte + p3_c[i]);
+        tte = qvp->eps_c[i];
+        bl1->eps[i] = g1*(qvp->eps_c[i] + b1*qvp->p3_c[i]);
+        bl1->p3[i] = g1*(b1*tte + qvp->p3_c[i]);
         sueps = sueps + bl1->eps[i];
       } //enddo
 
+      goto pnu987;
+      // With a potential for the projectile it is nonsens to search for 
+      // nucleons ON shell  31/05/2010.      
+
       cobe = (tlab + fmpinc)/sueps;
 
@@ -3012 +3472 @@
       if(iflag == nbtest) { // too much..all momentum to 0
         for(G4int klm = 1; klm <= bl3->ia1; klm++) { //do klm=1,bl3->ia1
-          eps_c[klm] = fmp;
+          qvp->eps_c[klm] = fmp;
           bl1->p1[klm] = 0.0;
           bl1->p2[klm] = 0.0;
-          p3_c[klm] = 0;
+          qvp->p3_c[klm] = 0;
         }
         goto pnu1870;
@@ -3030 +3490 @@
             }
           }
-          eps_c[i_emax] = fmp;
+          qvp->eps_c[i_emax] = fmp;
           bl1->p1[i_emax] = 0.0;
           bl1->p2[i_emax] = 0.0;
-          p3_c[i_emax] = 0.0;
+          qvp->p3_c[i_emax] = 0.0;
 
           if(i_emax == bl3->ia1) { //    circular permut if the last one
-            epsv = eps_c[bl3->ia1]; //   permutation circulaire
+            epsv = qvp->eps_c[bl3->ia1]; //      permutation circulaire
             p1v = bl1->p1[bl3->ia1];
             p2v = bl1->p2[bl3->ia1];
-            p3v = p3_c[bl3->ia1];
+            p3v = qvp->p3_c[bl3->ia1];
             for(bl2->k = bl3->ia1-1; bl2->k >= 1; bl2->k = bl2->k - 1) { //do k=bl3->ia1-1,1,-1
-              eps_c[bl2->k+1] = eps_c[bl2->k];
+              qvp->eps_c[bl2->k+1] = qvp->eps_c[bl2->k];
               bl1->p1[bl2->k+1] = bl1->p1[bl2->k];
               bl1->p2[bl2->k+1] = bl1->p2[bl2->k];
-              p3_c[bl2->k+1] = p3_c[bl2->k];
+              qvp->p3_c[bl2->k+1] = qvp->p3_c[bl2->k];
             }
-            eps_c[1] = epsv;
+            qvp->eps_c[1] = epsv;
             bl1->p1[1] = p1v;
             bl1->p2[1] = p2v;
-            p3_c[1] = p3v;      // fin permut.
+            qvp->p3_c[1] = p3v;         // fin permut.
           }
           sp1t1 = 0.0;   // re-compute the last one 
@@ -3057 +3517 @@
             sp1t1 = sp1t1 + bl1->p1[j];
             sp2t1 = sp2t1 + bl1->p2[j];
-            sp3t1 = sp3t1 + p3_c[j];
+            sp3t1 = sp3t1 + qvp->p3_c[j];
           }
           bl1->p1[bl3->ia1] = -sp1t1;
           bl1->p2[bl3->ia1] = -sp2t1;
-          p3_c[bl3->ia1] = -sp3t1;
-          eps_c[bl3->ia1] = w(bl1->p1[bl3->ia1],bl1->p2[bl3->ia1],p3_c[bl3->ia1],fmp);  
+          qvp->p3_c[bl3->ia1] = -sp3t1;
+          qvp->eps_c[bl3->ia1] = w(bl1->p1[bl3->ia1],bl1->p2[bl3->ia1],qvp->p3_c[bl3->ia1],fmp);        
 
           goto pnu1870;  // ..and boost all of them.
@@ -3083 +3543 @@
       }
       bl1->eps[ilm] = bl1->eps[ilm] + v0;  
-    }
-  }
+
+ pnu987:
+      //  G4cout <<"(checkpoint 'pnu987)" << G4endl;
+  // ici sauvegarde des energies et impulsions des nucleons du projectile i.l.:
+  // (hors potentiel)
+  for(G4int i = 1; i <= bl3->ia1; i++) {
+    qvp->eps_c[i]=bl1->eps[i];
+    qvp->p1_s[i]=bl1->p1[i];
+    qvp->p2_s[i]=bl1->p2[i];
+    qvp->p3_s[i]=bl1->p3[i];
+  }
+
+  // first particle entering in the target potential:   
+  bl1->eps[ilm]=bl1->eps[ilm]+v0;
+          
+  // correction en test 16/04/2009 (pourquoi pas la correction aussi de l'impulsion
+  // comme pour un nucleon in ???? ....on essaye de traiter les nucleons hors couche?)
+  // 29/4/2009: je pense que c'est une erreur et je transforme aussi maintenant l'impulsion:
+  var_ab = std::pow(bl1->p1[ilm],2)
+    + std::pow(bl1->p2[ilm],2)
+    + std::pow(bl1->p3[ilm],2);
+  if(var_ab > 0.) {
+    gpsg=std::sqrt((std::pow(bl1->eps[ilm],2)
+                    -fmp*fmp)/var_ab);
+  }
+  bl1->p1[ilm]=bl1->p1[ilm]*gpsg;             
+  bl1->p2[ilm]=bl1->p2[ilm]*gpsg;
+  bl1->p3[ilm]=bl1->p3[ilm]*gpsg;        
+    } // For pion-in
+  } // For nucleon-in
 
  pnu1871:
+  //  G4cout <<"(checkpoint 'pnu1871)" << G4endl;
+  //  bl1->dump(27);
+  //  bl3->dump();
+  //  randomGenerator->printSeeds();
   // evaluation of the times t(a,b)
   bl2->k = 0;
   kcol = 0;
-  if (kindstruct->kindf7 <= 2) {
+  if (kindstruct->kindf7 <= 2 && kindstruct->kindf7 > 0) {
     // modif s.vuillier tient compte propagation projectile,1e collision
     // imposee pour lui (c'est une maniere de faire!!)
@@ -3155 +3647 @@
   else {
     // deutons
-    if (kindstruct->kindf7 < 6) { //then
+    if (kindstruct->kindf7 < 6 && kindstruct->kindf7 > 0) { //then
       // here for incoming pions:
       for(G4int i = bl3->ia1+1; i <= ia; i++) { //do i=ia1+1,ia
@@ -3176 +3668 @@
     }
     else {
+      // Counting the spectators and transparents of composit projectiles:      
+      //      G4cout <<"(checkpoint 'setipszero1)" << G4endl;
+      ips=0;
+      nb_transprojo=0;
       for(G4int i = 1; i <= bl3->ia1; i++) { //do 38 i=1,ia1
         bl5->nesc[i] = 1;
+        // Spectators of composite projectiles (7/2006, AB)
+        // Pour garder trace des spectateurs du projectile:
+        npproj[i]=1;
+
         if (i != ilm) { 
           goto pnu36;
@@ -3193 +3693 @@
         goto pnu37;
       pnu36:
+        //      G4cout <<"(checkpoint 'pnu36)" << G4endl;
         t1 = bl3->x1[i]*(bl1->p1[i])+bl3->x2[i]*(bl1->p2[i])+bl3->x3[i]*(bl1->p3[i]);                      
         t2 = bl1->p1[i]*(bl1->p1[i])+bl1->p2[i]*(bl1->p2[i])+bl1->p3[i]*(bl1->p3[i]);               
@@ -3199 +3700 @@
         //   1379       c incoming nucleons enter potential at maximum radius (modif. 13/06/01)
         t5 = t3*t3 + ((ws->rmaxws)*(ws->rmaxws) - t4)/t2;
+        if (t5 < 0.) {
+          // this is a projectile spectator:
+          //      G4cout <<"(checkpoint 'incripspnu36)" << G4endl;
+          ips=ips+1;
+          isPartOfSpectatorNucleus[i] = true;
+          ps->n_projspec[ips] = i;
+          continue;
+          //goto pnu38;
+        } // endif
+                
         if(verboseLevel > 3) {
           G4cout <<"x1 = " << bl3->x1[i] <<" x2 = " << bl3->x2[i] <<" x3 = " << bl3->x3[i] << G4endl; 
@@ -3208 +3719 @@
           G4cout <<"t5 = " << t5 << G4endl;
         }
-        if (t5 < 0.) {
-          continue;
-        }
+        //      if (t5 < 0.) {
+        //        ips++;
+        //        ps->n_projspec[ips] = i;
+        //        continue;
+        //      }
         tref = (-1.0*t3 - std::sqrt(t5))*(bl1->eps[i]);  
-        if (tref > bl4->tmax5) {
-          continue;
-        }
+
+        // Drop reflection avatars with too large tref and the ones
+        // with negative tref (spectator, outside the nucleus)
+        //      if (tref > bl4->tmax5 || tref <= 0.0) {
+        //        continue;
+        //      }
+        if (tref > bl4->tmax5) { //then
+          // This is also a projectile spectator:
+          ips=ips+1;
+          ps->n_projspec[ips] = i;
+          continue; // go to 38
+        } // endif
+
+        // Test to avoid a nucleon entering in the past (3/06/2010 as in incl4.6 a.b.)
+        if(tref <= 0.0) { // then
+          // This is also a projectile spectator:
+          ips=ips+1;
+          ps->n_projspec[ips] = i;
+          continue; // go to 38
+        } // endif
+      
         npproj[i] = 1;
       pnu37:
+        //      G4cout <<"(checkpoint 'pnu37)" << G4endl;
         bl2->k = bl2->k + 1; 
         bl2->crois[bl2->k] = tref; 
         bl2->ind[bl2->k] = i; 
         bl2->jnd[bl2->k] = -1;
-      }
+      //    pnu38:
+      }
+
+      // For composit beams, cannot decide that there is no interaction,
+      // (another nucleon can enter later at RMAX).
       kcol = 1;
-      
+      n_activnuc = ia-ips;  //number of active nucleons (not missing the target)
+      n_enter_pot = ia-ips; //number of projectile nucleons that can enter
+
+      /*
+      if(useProjSpect) { // Treat the projectile spectators:
+        n_activnuc = ia - ips; // number of active nucleons (not missing the target)
+        if(ips != 0) { //then
+          e_spec = 0.0;
+          p1_spec = 0.0;
+          p2_spec = 0.0;
+          p3_spec = 0.0;
+          z_projspec = 0;
+          for(G4int i_spec = 1; i_spec <= ips; i_spec++) { // do i_spec=1,ips
+            G4int i_c = n_projspec[i_spec];
+            e_spec = e_spec + bl1->eps[i_c];
+            p1_spec = p1_spec + bl1->p1[i_c];
+            p2_spec = p2_spec + bl1->p2[i_c];
+            p3_spec = p3_spec + bl1->p3[i_c];
+            if(bl1->ind2[i_c] == 1) z_projspec = z_projspec + 1;
+          }
+          G4double p_spec2 = std::pow(p1_spec,2) + std::pow(p2_spec,2) + std::pow(p3_spec,2);
+          G4double s_spec = std::sqrt(std::pow(e_spec,2) - p_spec2);
+          // write(6,*) 'e,p,sqs:',e_spec,sqrt(p_spec2),s_spec
+
+          // c masses from dresner code:
+          G4int fja = ips;                                                                   
+          G4int fjz = z_projspec;                                                                   
+          G4double ex2 = abla->eflmac(fja, fjz, 0, 0); // dnrgy(fja,fjz);
+          G4double rnmass = ex2+um*fja;
+
+          // c write(6,*) fja,fjz,ex2,um                                                         
+          ex_projspec = s_spec - rnmass;
+          if(ex_projspec <= 0.0) { //then
+            // c no spectator nucleus compatible with available energy:
+            ips = 0;
+            // c n_activnuc = ia-ips!number of active nucleons (not missing the target)
+            goto pnu239;
+          }
+          t_projspec = e_spec - s_spec;
+          G4double coef = std::sqrt(t_projspec*(t_projspec + 2.0*s_spec)/p_spec2);
+          p1_projspec = p1_spec*coef;
+          p2_projspec = p2_spec*coef;
+          p3_projspec = p3_spec*coef;
+          m_projspec = rnmass;
+
+          // c check module...      
+          //    c      write(6,*) 'projectile spectator a,z:',a_projspec,z_projspec
+          //    c      write(6,*) 'm,t,e*:',m_projspec,t_projspec,ex_projspec
+          //    c      write(6,*) 'p1,p2,p3:',p1_projspec,p2_projspec,p3_projspec
+          //    c      write(6,*) '   constituants:',ips
+          //    c      do i_spec=1,ips
+          //    c        i_c = n_projspec(i_spec)
+          //    c      write(6,*) i_c,eps(i_c),p1(i_c),p2(i_c),p3(i_c)       
+          // c      enddo
+          // c ...end check module 
+        } // endif
+      pnu239:
+        a_projspec = ips; // number of nucleons in the spectator nucleus
+      } // End of projectile spectator treatment
+      */
+
       //      for(G4int i = bl3->ia1+1; i < ia; i++) { //do  39 i=ia1+1,ia
       for(G4int i = bl3->ia1+1; i <= ia; i++) { //do  39 i=ia1+1,ia
@@ -3298 +3894 @@
     goto pnu48;
   }
+  // Here transparent event (no interaction avatar found)
+  //  G4cout <<"(checkpoint 'heretransparent)" << G4endl;
   nopart = -1;
   // Pour eviter renvoi des resultats du run precedent cv 7/7/98
@@ -3316 +3914 @@
   // Initialization  at the beginning of the run
  pnu48:
+  //  G4cout <<"(checkpoint 'pnu48)" << G4endl;
   if(verboseLevel > 3) {
     G4cout <<"Beginning a run..." << G4endl;
@@ -3364 +3963 @@
   }
  pnu449:
+  //  G4cout <<"(checkpoint 'pnu449)" << G4endl;
   if(verboseLevel > 3) {
     G4cout <<"Now at 449" << G4endl;
@@ -3372 +3972 @@
 
  pnu44:
+  //  G4cout <<"(checkpoint 'pnu44)" << G4endl;
   if(verboseLevel > 3) {
     G4cout <<"Now at 44" << G4endl;
@@ -3406 +4007 @@
   }
  pnu448:
+  //  G4cout <<"(checkpoint 'pnu448)" << G4endl;
+  //  bl2->dump();
   imin = indic[next];
   bl9->l1 = bl2->ind[imin]; //NOTE: l1 changed to bl9->l1.
@@ -3419 +4022 @@
   // test le 20/3/2003: tue sinon le dernier avatar?
   if (bl2->k == 0) {
-    if(verboseLevel > 2) {
+    if(verboseLevel > -2) {
       G4cout <<"k == 0. Going to the end of the avatar." << G4endl;
     }
@@ -3438 +4041 @@
   }
  pnu46:
+  //  G4cout <<"(checkpoint 'pnu46)" << G4endl;
   if(verboseLevel > 3) {
     G4cout <<"G4Incl: Now at pnu46." << G4endl;
   }
   tim = timi + tau;
+  avatarCounter++;
+  //  G4cout <<"FOUND AVATAR " <<  avatarCounter << " : tau = " << tau;
+  //  G4cout <<"  l1 = " << bl9->l1 << "   l2 = " << bl9->l2 << G4endl;
 
   // tableau des energies a t=20,40,60 fm/c
@@ -3496 +4103 @@
     }
   }
+  //  G4cout <<"(checkpoint 'pnu645)" << G4endl;
+  // Here, stoping condition for heavy (light) ions:
+  if(kindstruct->kindf7 < 0) {
+    if(n_activnuc-nbquit <= ia_breakup) {
+      //      G4cout <<"stopped cascade for fermi breakup " << G4endl;
+      goto pnu255;
+    }
+  }
 
   // modif: pas de reflexions avant au moins un avatar du (des) nucleon incident
   // celui-ci ne peut etre qu'une collision nn (ou pin)
 
+  /*
   if((irst_avatar == 0) && (bl9->l2 == -1)) {
     if(verboseLevel > 3) {
@@ -3508 +4124 @@
 
   irst_avatar = irst_avatar+1;
+  */
+  //  if((bl3->ia1 + bl3->ia2 - nbquit) <= 10) goto pnu255; // Stop the cascade if only very few nucleons are left.
 
   if (tim < temfin) {
@@ -3518 +4136 @@
   goto pnu255;
  pnu49:
+  //  G4cout <<"(checkpoint 'pnu49)" << G4endl;
   if(verboseLevel > 3) {
     G4cout <<"G4Incl: Now at pnu49. " << G4endl;
@@ -3551 +4170 @@
   // l1 est un delta: 
  pnu830:
+  //  G4cout <<"(checkpoint 'pnu830)" << G4endl;
   if(verboseLevel > 3) {
     G4cout <<"G4Incl: Now at pnu830." << G4endl;
@@ -3566 +4186 @@
   }
  pnu803:
+  //  G4cout <<"(checkpoint 'pnu803)" << G4endl;
   // pas de collision entre 2 non participants:
   if(jparticip[bl9->l1] == 0 && jparticip[bl9->l2] == 0) {
@@ -3632 +4253 @@
   goto pnu261;
  pnu260:
+  //  G4cout <<"(checkpoint 'pnu260)" << G4endl;
   bmax2 = totalCrossSection(sq,mg,isos)/31.41592;
  pnu261:
+  //  G4cout <<"(checkpoint 'pnu261)" << G4endl;
   if (bb2 < bmax2) {
     goto pnu220;
@@ -3648 +4271 @@
   // evaluation of the positions at time = tim
  pnu220:
+  //  G4cout <<"(checkpoint 'pnu220)" << G4endl;
   timi = tim;
    if(verboseLevel > 3) {
@@ -3682 +4306 @@
 
   // gel des nucleons non participants sur le premier avatar (nn)=(l1,1)      
+  /*
   if (irst_avatar == 1) {
     for(G4int i = 1; i <= bl9->l1; i = i + bl9->l1 - 1) { // bugfix!
@@ -3697 +4322 @@
   }
   else {
+  */
     for(G4int i = 1; i <= ia; i++) {
       bl1->ta = tau/bl1->eps[i];
@@ -3703 +4329 @@
       bl3->x3[i] = bl3->x3[i] + bl1->p3[i]*(bl1->ta);
     }
-  }
+    //  }
 
   //  if(npion != 0) {
@@ -3714 +4340 @@
     }
   }
+  //  G4cout <<"(checkpoint 'pnu840)" << G4endl;
 
   if(bl9->l2 == 0) {
@@ -3828 +4455 @@
   }
  pnu243:
+  //  G4cout <<"(checkpoint 'pnu243)" << G4endl;
   if (bl1->ind1[bl9->l2] == 1) { // bugfix 1 -> 0
     goto pnu241;
@@ -3838 +4466 @@
   goto pnu241;
  pnu248:
+  //  G4cout <<"(checkpoint 'pnu248)" << G4endl;
   if(verboseLevel > 3) {
     G4cout <<"Pauli blocked transition!" << G4endl;
@@ -3874 +4503 @@
 
  pnu241:
-
+  //  G4cout <<"(checkpoint 'pnu241)" << G4endl;
   // la premiere collision a deux corps ne peut pas baisser l'energie
   // du nucleon de recul (bloque par pauli dans un noyau cible froid).
@@ -3974 +4603 @@
   q4[npion] = t[14]; //t(15)->t[14]
  pnu240:
+  //  G4cout <<"(checkpoint 'pnu240)" << G4endl;
   ncol = ncol + 1;
   if (bl9->l2 != 1) {
@@ -4045 +4675 @@
 
   pnu870:
+  //  G4cout <<"(checkpoint 'pnu870)" << G4endl;
   bl5->tlg[bl9->l1] = th*(bl1->eps[bl9->l1])/std::sqrt(std::pow(bl1->eps[bl9->l1],2)-std::pow(bl1->p1[bl9->l1],2)-std::pow(bl1->p2[bl9->l1],2)-std::pow(bl1->p3[bl9->l1],2));
   bl5->tlg[bl9->l2] = th*(bl1->eps[bl9->l2])/std::sqrt(std::pow(bl1->eps[bl9->l2],2)-std::pow(bl1->p1[bl9->l2],2)-std::pow(bl1->p2[bl9->l2],2)-std::pow(bl1->p3[bl9->l2],2));
@@ -4250 +4881 @@
 
  pnu848:
+  //  G4cout <<"(checkpoint 'pnu848)" << G4endl;
   if (npion == 0) {
     goto pnu844;
@@ -4269 +4901 @@
 
  pnu843:
+  //  G4cout <<"(checkpoint 'pnu843)" << G4endl;
   if(bl1->ind1[bl9->l2] != 1) {
     for(G4int k20 = 1; k20 <= npion; k20++) {
@@ -4276 +4909 @@
   }
  pnu844:
-
+  //  G4cout <<"(checkpoint 'pnu844)" << G4endl;
   if(bl1->ind1[bl9->l1]+bl1->ind1[bl9->l2] <= ich1+ich2) {
     goto pnu849;
@@ -4286 +4919 @@
   goto pnu821;
  pnu820:
+  //  G4cout <<"(checkpoint 'pnu820)" << G4endl;
   if(bl1->ind1[bl9->l2]-ich2 != 1) {
     goto pnu849;
@@ -4292 +4926 @@
 
  pnu821:
+  //  G4cout <<"(checkpoint 'pnu821)" << G4endl;
   standardRandom(&rndm,&(hazard->igraine[16]));
   // largeur variable du delta (phase space factor G4introduced 4/2001)
@@ -4320 +4955 @@
   // decay of the delta particle                                       p-n09780
  pnu805:
+  //  G4cout <<"(checkpoint 'pnu805)" << G4endl;
   npion = npion + 1;
   ichd = bl1->ind2[bl9->l1];
@@ -4372 +5008 @@
 
  pnu837:
+  //  G4cout <<"(checkpoint 'pnu837)" << G4endl;
   ipi[npion]=bl1->ind2[bl9->l1]*2;
   bl1->ind2[bl9->l1]=-1*(bl1->ind2[bl9->l1]);
   goto pnu809;
  pnu806:
+  //  G4cout <<"(checkpoint 'pnu806)" << G4endl;
   bl1->ind2[bl9->l1]=bl1->ind2[bl9->l1]/3;
   ipi[npion]=2*(bl1->ind2[bl9->l1]);
  pnu809: // continue
+  //  G4cout <<"(checkpoint 'pnu809)" << G4endl;
   bl1->ind1[bl9->l1]=0;
   bl5->tlg[bl9->l1]=0.;
@@ -4422 +5061 @@
   // it is blocked!      
  pnu1848:
+  //  G4cout <<"(checkpoint 'pnu1848)" << G4endl;
   mpaul2 = mpaul2 + 1;
   if(varavat->kveux == 1) {
@@ -4467 +5107 @@
   // valid decay of the delta
  pnu850:
+  //  G4cout <<"(checkpoint 'pnu850)" << G4endl;
   if(varavat->kveux == 1) {
     varavat->bloc_paul[iavat] = 0;
@@ -4514 +5155 @@
   
  pnu4047:
+  //  G4cout <<"(checkpoint 'pnu4047)" << G4endl;
   if (bl5->nesc[bl9->l1] != 0) {
     goto pnu845;
@@ -4538 +5180 @@
 
  pnu845: 
+  //  G4cout <<"(checkpoint 'pnu845)" << G4endl;
   new2(y1[npion], y2[npion], y3[npion], q1[npion], q2[npion], q3[npion], q4[npion], npion, bl9->l1);
   if(bl2->k == 0) {
@@ -4550 +5193 @@
   // pion-nucleon collision
  pnu801:
+  //  G4cout <<"(checkpoint 'pnu801)" << G4endl;
   if(verboseLevel > 3) {
     G4cout <<"Pion-nucleon collision!" << G4endl;
@@ -4588 +5232 @@
 
  pnu832:
+  //  G4cout <<"(checkpoint 'pnu832)" << G4endl;
   if (bl2->k == 0) {
     goto pnu230;
@@ -4597 +5242 @@
   goto pnu44;
  pnu831:
+  //  G4cout <<"(checkpoint 'pnu831)" << G4endl;
   standardRandom(&rndm, &(hazard->igraine[18]));
   geff = t[9]/sq; //t(10)->t[9]
@@ -4714 +5360 @@
   // reflection on or transmission through the potential wall
  pnu600:
+  //  G4cout <<"(checkpoint 'pnu600)" << G4endl;
   // deutons pas bien compris ici cv ?
   if (npproj[bl9->l1] == 0) {
@@ -4736 +5383 @@
   }
   
+  if (bl1->ind2[bl9->l1] == 1) { //then
+    npenter=npenter+1;
+  } else {
+    nnenter=nnenter+1;
+  } // endif
+
   var_ab = std::pow(bl1->p1[bl9->l1],2) + std::pow(bl1->p2[bl9->l1],2) + std::pow(bl1->p3[bl9->l1],2);
   gpsg = 0.0;
@@ -4756 +5409 @@
   // pour un non participant la transmission est impossible:
  pnu608:
+  //  G4cout <<"(checkpoint 'pnu608)" << G4endl;
   if(varavat->kveux == 1) {
     varavat->del1avat[iavat] = bl1->ind1[bl9->l1];
@@ -4777 +5431 @@
 
  pnu605:
+  //  G4cout <<"(checkpoint 'pnu605)" << G4endl;
   fm = fmp;
   pot = v0;
 
  pnu606:
+  //  G4cout <<"(checkpoint 'pnu606)" << G4endl;
   if(verboseLevel > 3) {
     G4cout <<"G4Incl: Now at pnu606. Calculating transmission probability." << G4endl;
   }
-  tp = transmissionProb(bl1->eps[bl9->l1]-fm,bl1->ind2[bl9->l1],itch,bl3->r2,v0);
+  tp = transmissionProb(bl1->eps[bl9->l1]-fm,bl1->ind2[bl9->l1],1,itch,bl3->r2,v0);
   if(varavat->kveux == 1) {
     varavat->energyavat[iavat] = bl1->eps[bl9->l1] - fm;
@@ -4796 +5452 @@
     goto pnu601;
   }
-
+  // Commande clusters ou non (1/0)! 
+  if(kclst == 0) goto pnu1610;
+  // Cluster goes here!
+ pnu1610:
+  //  G4cout <<"(checkpoint 'pnu1610)" << G4endl;
   // ici la particule l1 s'échappe du noyau:
   bl5->nesc[bl9->l1] = 1;
@@ -4811 +5471 @@
   bl1->eps[bl9->l1] = bl1->eps[bl9->l1] - pot;
 
+  // Test for a transparent nucleon for composit beams and stored for FermiBreakup
+  if(bl3->ia1 > 1 && nc[bl9->l1] == 0) {
+    //    G4cout <<"(checkpoint 'incripspnu1610)" << G4endl;
+    ips=ips+1;
+    isPartOfSpectatorNucleus[bl9->l1] = true;
+    ps->n_projspec[ips]=bl9->l1;
+    nb_transprojo = nb_transprojo + 1;
+  } // endif
+
+  // pnu610: goes here
+
   // comptage des particules hors du noyau (7/6/2002):
   // (remnant minimum=1 nucleon)
@@ -4827 +5498 @@
   // here no transmission possible
  pnu601:
+  //  G4cout <<"(checkpoint 'pnu601)" << G4endl;
   pspr=bl3->x1[bl9->l1]*(bl1->p1[bl9->l1])+bl3->x2[bl9->l1]*(bl1->p2[bl9->l1])+bl3->x3[bl9->l1]*(bl1->p3[bl9->l1]);
   if(varavat->kveux == 1) {
@@ -4841 +5513 @@
 
  pnu602:
+  //  G4cout <<"(checkpoint 'pnu602)" << G4endl;
   if(verboseLevel > 3) {
     G4cout <<"G4Incl: Now at pnu602." << G4endl;
   }
   
-  if(bl2->k != 0) {
+  if(bl2->k != 0) { // GOTO elimination?
     kd = 0;
     ccr = tau;
@@ -4862 +5535 @@
     bl2->k = bl2->k - kd;
 
-    if (bl5->nesc[bl9->l1] == 1) {
+    if (bl5->nesc[bl9->l1] != 0) { // modif 10/02 pour logique des clusters (nesc()=2!)
       goto pnu613;
     }
@@ -4898 +5571 @@
 
  pnu615:
+  //  G4cout <<"(checkpoint 'pnu615)" << G4endl;
   if(verboseLevel > 3) {
     G4cout <<"G4Incl: Now at pnu615." << G4endl;
@@ -4912 +5586 @@
   }
  pnu613:
+  //  G4cout <<"(checkpoint 'pnu613)" << G4endl;
   if(verboseLevel > 3) {
     G4cout <<"G4Incl: Now at pnu613." << G4endl;
@@ -4948 +5623 @@
   // decay of the surviving deltas
  pnu230:
+  //  G4cout <<"(checkpoint 'pnu230)" << G4endl;
   if(verboseLevel > 3) {
     G4cout <<"G4Incl: Now at pnu230." << G4endl;
@@ -4953 +5629 @@
   // decay of the surviving deltas
 pnu255:
+  //  G4cout <<"(checkpoint 'pnu255)" << G4endl;
   if(verboseLevel > 3) {
     G4cout <<"G4Incl: Now at pnu6255." << G4endl;
@@ -4965 +5642 @@
 
   npidir = npion;
+
+  idq = 0;
+  destar = 0.;
+
   for(G4int i = 1; i <= ia; i++) {
+    G4int iblcdpp=0;
     if (bl1->ind1[i] == 0) {
       continue;
@@ -4980 +5662 @@
     xy3 = bl1->p3[i];
     xye = bl1->eps[i];
+    ichd=bl1->ind2[i];
     if(varavat->kveux == 1) {
       iavat = iavat + 1;
@@ -5000 +5683 @@
       G4cout <<"i = " << i << " bl1->p1[i] = " << bl1->p1[i] << " bl1->p2[i] = " << bl1->p2[i] << " bl1->p3[i] = " << bl1->p3[i] << " bl1->eps[i] = " << bl1->eps[i] << G4endl;
     }
-    
+
+// C 6/08/09; Don't have to check CDPP for the forced decay of a DELTA OUTSIDE!
+//       IF(NESC(I).EQ.1) GO TO 1850
+      
+// C------
+// C On teste CDPP apres la decroissance forcee June/2005 (AB-CV)
+// C    Si rejet de la decroissance -> delta devient nucleon sans pion 
+// C       et toute l'energie en exces est mise en energie d'excitation du noyau
+// C       On respecte aussi la conservation de charge (Z du remnant!)
+ 
+// C Le decay ne peut conduire a un noyau de A nucleons
+// C  sous l'energie de Fermi et dans une config. d'energie inferieure a
+// C  EFER-(IA2-NBALTTF)*TF).
+//        EGS=0.
+//        NBALTTF=0
+//        DO II=1,IA
+//          IF(NESC(II).EQ.0) THEN
+//          IF(SQRT(P1(II)**2+P2(II)**2+P3(II)**2).LT.PF) THEN
+//             NBALTTF=NBALTTF+1
+//             EGS=EGS+EPS(II)-FMP
+//          ENDIF
+//       ENDIF
+//        END DO
+//        IF(EGS.GE.(EFER-(IA2-NBALTTF)*TF)) GO TO 1850
+// C ATTENTION, logique negative!!! Liberer le goto si on veut supprimer la
+// C   sequence precedente (CDPP sur Delta-> pi N)
+// C        GO TO 1850
+
+// C Ici ce decay est interdit par CDPP.
+// C Restitution du Delta: un neutron au repos et correction conservee pour 
+// C la charge IDQ et l'energie d'excitation DESTAR du remnant. On viole donc
+// C un peu la conservation d'impulsion, mais semble moins grave que de convertir
+// C la masse d'un pion en impulsion d'un nucleon unique....
+//       IBLCDPP=1              !flag de blocage par CDPP
+//       IND1(i) = 0    !identif delta (nucleon)
+//       IND2(i) = -1   !isospin du neutron
+//      IDQ = IDQ + (ICHD+1)/2
+//       p1(i)=0.               !impulsion energie
+//       p2(i)=0.
+//       p3(i)=0.
+//       eps(i)=FMP
+//      DESTAR = DESTAR + xye - FMP
+//       NPION=NPION-1  !pas le dernier pion
+
+// 1850  CONTINUE
+// C------    
     if(bl5->nesc[i] == 0) {
       idecf = 1;
@@ -5022 +5750 @@
       // fin surface
     }
+    if(iblcdpp == 1) continue; // goto pnu257;
 
     if (bl1->ind2[i]*(bl1->ind2[i]) == 9) {
@@ -5036 +5765 @@
 
   pnu283:
+    //    G4cout <<"(checkpoint 'pnu283)" << G4endl;
     if(verboseLevel > 3) {
       G4cout <<"G4Incl: Now at pnu283." << G4endl;
@@ -5045 +5775 @@
 
   pnu280:
+    //    G4cout <<"(checkpoint 'pnu285)" << G4endl;
     if(verboseLevel > 3) {
       G4cout <<"G4Incl: Now at pnu280." << G4endl;
@@ -5053 +5784 @@
 
   pnu285:
+    //    G4cout <<"(checkpoint 'pnu285)" << G4endl;
     if(verboseLevel > 3) {
       G4cout <<"G4Incl: Now at pnu285." << G4endl;
@@ -5094 +5826 @@
 
  pnu256:
+  //  G4cout <<"(checkpoint 'pnu256)" << G4endl;
   if(verboseLevel > 3) {
     G4cout <<"G4Incl: Now at pnu256." << G4endl;
@@ -5102 +5835 @@
   npx = 0;
   erem = 0.;
+  // Excitation energy of the final delta of blocked decay: (AB CV 6/2005)
+  erem = erem + destar;
   izrem = 0;
+  // Correct the charge of the remnant if the final delta is blocked: (A.B., C.V. 6/2005)
+  izrem = izrem + idq;
   inrem = 0;
   iarem = 0;
@@ -5115 +5852 @@
   pout3 = 0.0;
   eout = 0.0;
+  p1spec=0.;
+  p2spec=0.;
+  p3spec=0.;
   cmultn = 0.0;
 
-  if (kindstruct->kindf7 <= 2) {
+  if (kindstruct->kindf7 <= 2 && kindstruct->kindf7 > 0) {
     if (ncol == 0 || nc[1] == 0) { // then nc(1)->nc[0]
       if(verboseLevel > 3) {
@@ -5126 +5866 @@
   }
   else {
-    if (kindstruct->kindf7 <= 5) {
+    if (kindstruct->kindf7 <= 5 && kindstruct->kindf7 > 0) {
       if (ncol == 0) {
         if(verboseLevel > 3) {
@@ -5149 +5889 @@
   // pour eviter renvoi des resultats du run precedent cv 20/11/98
  pnu9100:
+  //  G4cout <<"(checkpoint 'pnu9100)" << G4endl;
   iarem = bl3->ia2;
   izrem = iz2;
@@ -5161 +5902 @@
 
  pnu9101:
+  //  G4cout <<"(checkpoint 'pnu9101)" << G4endl;
   if(verboseLevel > 3) {
     G4cout <<"G4Incl: Now at pnu9101." << G4endl;
@@ -5168 +5910 @@
   ekout = 0.0;
 
-  for(G4int i = 1; i <= ia; i++) {
-    if(bl5->nesc[i] != 0) {
-      xl1 = xl1-bl3->x2[i]*(bl1->p3[i]) + (bl3->x3[i])*(bl1->p2[i]);
-      xl2 = xl2-bl3->x3[i]*(bl1->p1[i]) + bl3->x1[i]*(bl1->p3[i]);
-      xl3 = xl3-bl3->x1[i]*(bl1->p2[i]) + bl3->x2[i]*(bl1->p1[i]);
-
-      // ici ajout de pout cv le 5/7/95
-      pout1 = pout1 + bl1->p1[i];
-      pout2 = pout2 + bl1->p2[i];
-      pout3 = pout3 + bl1->p3[i];
-      eout = eout+bl1->eps[i] - fmp;
-      //      ic33 = int((std::floor(double(bl1->ind2[i]+3))/double(2)));
-      ic33 = (bl1->ind2[i]+3)/2;
-      //nopart = nopart + 1;
-      kind[nopart] = 3 - ic33;
-      if(verboseLevel > 3) {
-        G4cout <<"kind[" << nopart << "] = " << kind[nopart] << G4endl;
-      }
-      ep[nopart] = bl1->eps[i] - fmp;
-      if(verboseLevel > 2) {
-        if(ep[nopart] > calincl->f[2]) {
-          G4cout <<"ep[" << nopart << "] = " << ep[nopart] << " > " << calincl->f[2] << G4endl;
-          G4cout <<"E = " << bl1->eps[nopart] << G4endl;
-          G4cout <<"px = " << bl1->p1[nopart] << " py = " << bl1->p2[nopart] << " pz = " << bl1->p3[nopart] << G4endl;
-          G4cout <<"Particle mass = " << fmp << G4endl;
-        }
-      }
-      bmass[nopart] = fmp;
-      ekout = ekout + ep[nopart];
-      ptotl = std::sqrt(std::pow(bl1->p1[i],2) + std::pow(bl1->p2[i],2) + std::pow(bl1->p3[i],2));
-      alpha[nopart] = bl1->p1[i]/ptotl;
-      beta[nopart] = bl1->p2[i]/ptotl;
-      gam[nopart] = bl1->p3[i]/ptotl;
-      nopart = nopart + 1;
-      continue;
-    }
-    
-    t[3] = bl3->x1[i]*(bl3->x1[i]) + bl3->x2[i]*(bl3->x2[i]) + bl3->x3[i]*(bl3->x3[i]); //t(4)->t[3]
-    erem = erem + bl1->eps[i] - fmp;
-    rcm1 = rcm1 + bl3->x1[i];
-    rcm2 = rcm2 + bl3->x2[i];
-    rcm3 = rcm3 + bl3->x3[i];
-    prem1 = prem1 + bl1->p1[i];
-    prem2 = prem2 + bl1->p2[i];
-    prem3 = prem3 + bl1->p3[i];
-    izrem = izrem + (1 + bl1->ind2[i])/2;
-    iarem = iarem + 1;
-  }
-
+  if(kclst != 0) { // Clusters go here!
+  }
+
+// Treat here the projectile spectators (including transparents):
+  projo_spec(bl3->ia1,ips,fmpinc,pinc,tlab);
+
+  for(G4int i = 1; i <= ia; i++) { // do 258 i=1,ia
+    //  pnu259:
+    if (bl5->nesc[i] == 0) goto pnu254;
+    if (bl5->nesc[i] > 1) continue; // on evite les nucleons des clusters
+
+    // here nesc=1 or negative value (emitted nucleons, spectators included)      
+    xl1_i=bl3->x2[i]*bl1->p3[i]-bl3->x3[i]*bl1->p2[i];  // moment angulaire emporte
+    xl2_i=bl3->x3[i]*bl1->p1[i]-bl3->x1[i]*bl1->p3[i];
+    xl3_i=bl3->x1[i]*bl1->p2[i]-bl3->x2[i]*bl1->p1[i];
+
+    xl1=xl1-xl1_i;                                                    
+    xl2=xl2-xl2_i;                                                     
+    xl3=xl3-xl3_i;                                                     
+
+    //   ici ajout de pout cv le 5/7/95
+    pout1=pout1+bl1->p1[i];
+    pout2=pout2+bl1->p2[i];
+    pout3=pout3+bl1->p3[i];
+    //      write(6,*)'eout nucleon',eout
+    eout = eout + bl1->eps[i]-fmp;
+
+    ic33=(bl1->ind2[i]+3)/2;
+
+    // here we don't copy the spectators of the projectile bound in a
+    //   spectator nucleus (treated by the deexcitation module):
+    if(kindstruct->kindf7 < 0 && i <= bl3->ia1 && ps->a_projspec != 0) {
+      for(G4int iloc =  1; iloc <= ps->a_projspec; iloc++) { //do iloc=1,a_projspec
+        if(i == ps->n_projspec[iloc]) {
+          p1spec=p1spec+bl1->p1[i];
+          p2spec=p2spec+bl1->p2[i];
+          p3spec=p3spec+bl1->p3[i];
+          // for balance, a transparent from projectile is out:
+          if(i > ips-nb_transprojo) ekout=ekout+bl1->eps[i]-fmp;
+          goto pnu258_end_of_outer_loop;
+          //      continue;
+        } // endif
+      } // enddo
+    } // endif
+    //    if (isPartOfSpectatorNucleus[i]) goto pnu258_end_of_outer_loop;
+
+    nopart=nopart+1;
+
+    ekout=ekout+bl1->eps[i]-fmp;
+
+    kind[nopart]=3-ic33;
+    // spectators of composite projectiles (7/2006, ab)
+    if(npproj[i] == 1) kind[nopart] = -kind[nopart];
+
+    ep[nopart]=bl1->eps[i]-fmp;
+    bmass[nopart]=fmp;
+    ptotl=std::sqrt(std::pow(bl1->p1[i],2) + std::pow(bl1->p2[i],2)+std::pow(bl1->p3[i],2));
+    alpha[nopart]=bl1->p1[i]/ptotl;
+    beta[nopart]=bl1->p2[i]/ptotl;
+    gam[nopart]=bl1->p3[i]/ptotl;
+ 
+    continue;
+  pnu254:
+    t[4]=bl3->x1[i]*bl3->x1[i]+bl3->x2[i]*bl3->x2[i]+bl3->x3[i]*bl3->x3[i];
+    erem=erem+bl1->eps[i]-fmp;
+    rcm1=rcm1+bl3->x1[i];
+    rcm2=rcm2+bl3->x2[i];
+    rcm3=rcm3+bl3->x3[i];
+    prem1=prem1+bl1->p1[i];
+    prem2=prem2+bl1->p2[i];
+    prem3=prem3+bl1->p3[i];
+    izrem=izrem+(1+bl1->ind2[i])/2;
+    iarem=iarem+1;
+  pnu258_end_of_outer_loop:
+    continue;
+  }
+  //  258 continue  
+
+  //pnu258:
+  //  G4cout <<"(checkpoint 'pnu258)" << G4endl;
   //  correction pions 21/3/95 jc
   ichpion = 0;
@@ -5229 +6002 @@
       xl3 = xl3 - y1[ipion]*q2[ipion] + y2[ipion]*q1[ipion];
       ichpion = ichpion + ipi[ipion]/2;
-      //      nopart = nopart + 1;
+      nopart = nopart + 1;
       kind[nopart] = 4 - ipi[ipion]/2;
       ptotl = std::sqrt(std::pow(q1[ipion],2) + std::pow(q2[ipion],2) + std::pow(q3[ipion],2));
@@ -5238 +6011 @@
       beta[nopart] = q2[ipion]/ptotl;
       gam[nopart] = q3[ipion]/ptotl;
-      nopart++;
+      //      nopart++;
     }
   }
@@ -5354 +6127 @@
     }
 
-    for(G4int ipart = 0; ipart < nopart; ipart++) {
+    for(G4int ipart = 0; ipart <= nopart; ipart++) {
       ep[ipart] = ep[ipart]*fffc;
     }
@@ -5361 +6134 @@
   // modif boudard juillet 99 (il faut tenir compte de la renormalisation
   // des energies pour les impulsions.)
-  pfrem1 = 0.0;
-  pfrem2 = 0.0;
-  pfrem3 = pinc;
-  for(G4int ipart = 0; ipart < nopart; ipart++) {
-    xmodp = std::sqrt(ep[ipart]*(2.0*bmass[ipart] + ep[ipart]));
+  pfrem1= -p1spec;
+  pfrem2= -p2spec;
+  pfrem3=pinc-p3spec;
+//   if(useProjSpect) {
+//     pfrem1 = -p1_spec;
+//     pfrem2 = -p2_spec;
+//     pfrem3 = pinc - p3_spec;
+//   } else {
+//     pfrem1 = 0.0;
+//     pfrem2 = 0.0;
+//     pfrem3 = pinc;
+//   }
+  for(G4int ipart = 0; ipart <= nopart; ipart++) {
+    G4double xmodp = 0.0;
+    if(ep[ipart] < 0.0) { // Safeguard against particles with negative energy
+      xmodp = 0.0;
+      G4cout <<"The energy of particle " << ipart << " is negative (" << ep[ipart] << "). Forcing xmodp = " << xmodp << "." << G4endl;
+    } else {
+      xmodp = std::sqrt(ep[ipart]*(2.0*bmass[ipart] + ep[ipart]));
+    }
     pfrem1 = pfrem1 - alpha[ipart]*xmodp;
     pfrem2 = pfrem2 - beta[ipart]*xmodp;
@@ -5394 +6182 @@
     esrem = 0.0;
   }
+  goto pnureturn;
 
   if(verboseLevel > 3) {
@@ -5404 +6193 @@
   }
  pnureturn:
+  //  G4cout <<"(checkpoint 'pnureturn)" << G4endl;
+//   if(useProjSpect) {
+//     varntp->spectatorA = a_projspec;
+//     varntp->spectatorZ = z_projspec;
+//     varntp->spectatorEx = ex_projspec;
+//     varntp->spectatorT = t_projspec;
+//     varntp->spectatorM = m_projspec;
+//     varntp->spectatorP1 = p1_projspec;
+//     varntp->spectatorP2 = p1_projspec;
+//     varntp->spectatorP3 = p1_projspec;
+//   }
   (*ibert_p) = ibert;
   (*nopart_p) = nopart;
@@ -5415 +6215 @@
   (*bimpact_p) = bimpact;
   (*l_p) = l;
-}
-
+  (*xjrem) = 0;
+  (*yjrem) = 0;
+  (*zjrem) = 0;
+}
 
 void G4Incl::collis(G4double *p1_p, G4double *p2_p, G4double *p3_p, G4double *e1_p, G4double *pout11_p, G4double *pout12_p, 
@@ -6567 +7369 @@
     // pauli strict
     pmod = std::sqrt(std::pow(bl1->p1[l],2) + std::pow(bl1->p2[l],2) + std::pow(bl1->p3[l],2));
-    if (pmod < 270.0) {
+    if (pmod < fermi->pf) {
       return 1.0;
     } else {
@@ -6674 +7476 @@
   G4double sine = 0.0;
 
-  if((m-1) < 0) {
+  if((m-1) < 0) { // Nucleon - Nucleon
     sine = deltaProductionCrossSection(E,int(i));
   }
 
-  if((m-1) == 0) {
+  if((m-1) == 0) { // Nucleon - Delta
     sine = srec(E,(bl6->xx10),i,int(bl6->isa));
   }
 
-  if((m-1) > 0) {
+  if((m-1) > 0) { // Delta - Delta
     sine = 0.0;
   }
 
   stotResult = sine + lowEnergy(E,m,i);
+#ifdef G4INCLDEBUG
+  theLogger->fillHistogram1D("totalCXResult", stotResult);
+  theLogger->fillHistogram1D("totalLowECX", lowEnergy(E, m, i));
+  if(m == 0) { // Nucleon - Nucleon
+    if(i == 2) { // pp
+      theLogger->fillHistogram2D("ppCX", E, stotResult);
+    } else if(i == 0) { // pn
+      theLogger->fillHistogram2D("pnCX", E, stotResult);
+    } else if(i == -2) { // nn
+      theLogger->fillHistogram2D("nnCX", E, stotResult);
+    }
+  } else if(m == 1) { // Nucleon - Delta
+    theLogger->fillHistogram2D("NDCX", E, stotResult);
+  } else if(m == 2) { // Delta - Delta
+    theLogger->fillHistogram2D("DDCX", E, stotResult);
+  }
+#endif
   return stotResult;
 }
@@ -6858 +7677 @@
 }
 
-G4double G4Incl::transmissionProb(G4double E, G4double iz, G4double izn, G4double r, G4double v0)
-{
-  // transmission probability for a nucleon of kinetic energy
-  // E on the edge of the well of depth v0 (nr approximation)
-  // iz is the isospin of the nucleon,izn the instanteneous charge
-  // of the nucleus and r is the target radius
-
-  G4double x = 0.0;
+//      FUNCTION BARR(E,IZ,IA,IZN,R,V0)
+//  tp =          transmissionProb(bl1->eps[bl9->l1]-fm,bl1->ind2[bl9->l1],1,itch,bl3->r2,v0);
+G4double G4Incl::transmissionProb(G4double E, G4int iz, G4int ia, G4int izn,G4double R, G4double v0)
+{
+// C
+// C
+// C       NOUVEAU BARR RELATIVISTE ATTENTION AUX ENTREES SORTIES!
+// C    (06/2005)
+// C                                                                       P-N24450
+// CCC   BARR=TRANSMISSION PROBABILITY FOR A PARTICLE (Nucleon, cluster or 
+// CCC   pion) of kinetic energy E on the edge of the (attractive) well of 
+// CCC   depth V0 felt by the particle.
+// CCC          IZ is the isospin of the particle,
+// CCC          IZN the instantaneous charge of the nucleus and R the radius of
+// CCC       the well.
+// CCC          IA is the mass number of the particle, and MRE its mass energy.             
+// C                                                                       P-N24500
+// C     Modified 9/10/2002 for clusters (d,t,3He,4He) (IZ=isospin,IA=A)
+// C     IZ must be correct so that charge Q=(IA+IZ)/2.
+// C     Modified 4/2004 for relativistic expressions and pions.
+
+  G4double mre;
   G4double barr = 0.0;
-
-  // We need enough energy to escape from the potential well.
-  if (E > v0) {
-    x = std::sqrt(E*(E - v0));
-    barr = 4.0*x/(E + E - v0 + x + x);
-    if (iz > 0 && izn != 0) { // izn = 0 causes division by zero
-      G4double b = izn*1.44/r;
-      G4double px = std::sqrt((E - v0)/b);
-      
-      if (px < 1.0) {
-        G4double g = izn/137.03*std::sqrt(2.0*938.3/(E - v0))*(std::acos(px) - px*std::sqrt(1.0 - px*px));
-        if (g > 35.){
-          barr=0.0;
-        } else {
-          barr = barr*std::exp(-2.0*g);
-        }
-        return barr;
-      } else {
-        return barr;
-      }
-    } else {
-      return barr;
-    }
+  G4double v0ia = v0*ia;
+  mre = ia*938.3;
+
+  if (E > v0ia && izn > 0 && R != 0) goto tp1;
+  barr = 0.0;
+  return barr;
+
+ tp1:
+  G4double x=std::sqrt((2.*mre*E+E*E)*(2.*mre*(E-v0ia)+std::pow((E-v0ia),2)));                                          
+  barr=4.*x/(2.*mre*(2.*E-v0ia)+std::pow(E,2)+std::pow(E-v0ia,2)+2.*x);
+  G4double iq=(ia+iz)/2.;
+  if (iq > 0 && E > v0ia) goto tp2;
+  return barr;
+
+ tp2:
+  G4double b=iq*izn*1.44/R;                                             
+  G4double px=std::sqrt((E-v0ia)/b);                                             
+  if (px < 1.) goto tp3;
+  return barr;
+
+ tp3:
+  G4double g=iq*izn/137.03*std::sqrt(2.*mre/(E-v0ia)/(1.+(E-v0ia)/2./mre))
+    *(std::acos(px)-px*std::sqrt(1.-px*px));
+  if (g > 35.) {
+    barr=0.0;
   } else {
-    return barr;
-  }
+    barr=barr*std::exp(-2.*g);
+  } // endif
+  return barr;
 }
 
@@ -6908 +7744 @@
   G4double s_l = 0.0;
   G4double t1 = 0.0, t3 = 0.0, t4 = 0.0, t5 = 0.0;
-  if (ws->nosurf <= 0) { // modif pour w.s.:
-    xv = p/pf;
-    if(t2 <= pf2) {
-      r = interpolateFunction(xv);
-    } else {
-      r = ws->rmaxws;
-    }
-    r = r*r;
-  }
+   if (ws->nosurf <= 0) { // modif pour w.s.:
+     xv = p/pf;
+     if(t2 <= pf2) {
+       r = interpolateFunction(xv);
+#ifdef G4INCLDEBUG
+     theLogger->fillHistogram1D("refInterpolXV", xv);
+     theLogger->fillHistogram1D("refInterpolRes", r);
+     if(xv > 1.0) theLogger->fillHistogram1D("refInterpolResWhenXVgt1", r);
+#endif
+     } else {
+       r = ws->rmaxws;
+     }
+     r = r*r;
+   }
+/*
+  if (ws->nosurf <= 0) {
+    // modif pour w.s.:
+    xv=p/pf;
+    r=interpolateFunction(xv);
+    r=r*r;
+    if (t2 > pf2) r=std::pow(ws->rmaxws,2);
+  } //endif
+*/
  ref21:
   t4 = x1*x1 + x2*x2 + x3*x3;
@@ -6969 +7819 @@
   }
 
-  bl3->ia2 = int(std::floor(calincl->f[0] + 0.1)); // f(1) -> f[0]
+  //  bl3->ia2 = int(std::floor(calincl->targetA() + 0.1)); // f(1) -> f[0]
+  bl3->ia2 = calincl->targetA();
   sep = 6.8309;
 
@@ -6976 +7827 @@
       itg = 6 - 1;
     }
-    if(bl3->ia2 == 3 && calincl->f[1] == 1) {
+    if(bl3->ia2 == 3 && calincl->targetZ() == 1) {
       itg = 7 - 1;
     }
-    if(bl3->ia2 == 3 && calincl->f[1] == 2) {
+    if(bl3->ia2 == 3 && calincl->targetZ() == 2) {
       itg = 8 - 1;
     }
@@ -6988 +7839 @@
   }
 
-  if((calincl->f[2] >= 10.0) && (calincl->f[2] <= 100.0)) {
-    if(calincl->f[6] == 1.0) {
+  if((calincl->bulletE() >= 10.0) && (calincl->bulletE() <= 100.0)) {
+    if(calincl->bulletType() == 1) {
       bl3->ia1 = int(1.0);
       iz1 = 1.0;
       G4double fmpinc = 938.2796;
-      G4double pbeam2 = calincl->f[2]*(calincl->f[2] + 2.0*fmpinc);
+      G4double pbeam2 = calincl->bulletE()*(calincl->bulletE() + 2.0*fmpinc);
       bmaxt = ws->bmax;
-      proba_trans = coulombTransm(calincl->f[2],bl3->ia1,iz1,calincl->f[0],calincl->f[1]);
-      proba = forceAbs(1,calincl->f[0],calincl->f[1],calincl->f[2],bmaxt,proba_trans);
+      proba_trans = coulombTransm(calincl->bulletE(),bl3->ia1,iz1,calincl->targetA(),calincl->targetZ());
+      proba = forceAbs(1,calincl->targetA(),calincl->targetZ(),calincl->bulletE(),bmaxt,proba_trans);
       
       standardRandom(&alea,&(hazard->igraine[4]));
@@ -7003 +7854 @@
       }
 
-      (*iarem) = int(calincl->f[0]) + bl3->ia1;
-      (*izrem) = int(calincl->f[1]) + int(iz1);
+      (*iarem) = int(calincl->targetA()) + bl3->ia1;
+      (*izrem) = int(calincl->targetZ()) + int(iz1);
       
-      del = std::sqrt(std::pow(((calincl->f[0] + 1.0)*fmpinc + calincl->f[2]),2) - pbeam2);
+      del = std::sqrt(std::pow(((calincl->targetA() + 1.0)*fmpinc + calincl->bulletE()),2) - pbeam2);
       
-      (*erecrem) = pbeam2/((calincl->f[0] + 1.0)*fmpinc+calincl->f[2] + del);
-
-      (*esrem) = calincl->f[2] + sep - (*erecrem);
+      (*erecrem) = pbeam2/((calincl->targetA() + 1.0)*fmpinc+calincl->bulletE() + del);
+
+      (*esrem) = calincl->bulletE() + sep - (*erecrem);
 
       (*alrem) = 0.00001;
@@ -7019 +7870 @@
       return;
     }
-    else if((calincl->f[6] == 2) && (calincl->f[2] >= 20.0)) {
+    else if((calincl->bulletType() == 2) && (calincl->bulletE() >= 20.0)) {
       bl3->ia1 = int(1.0);
       iz1 = 0.0;
       G4double fmpinc = 938.2796;
-      G4double pbeam2 = calincl->f[2]*(calincl->f[2] + 2.0*fmpinc);
+      G4double pbeam2 = calincl->bulletE()*(calincl->bulletE() + 2.0*fmpinc);
       bmaxt = ws->bmax;
-      proba_trans = coulombTransm(calincl->f[2],bl3->ia1,iz1,calincl->f[0],calincl->f[1]);
-      proba = forceAbs(1,calincl->f[0],calincl->f[1],calincl->f[2],bmaxt,proba_trans);
+      proba_trans = coulombTransm(calincl->bulletE(),bl3->ia1,iz1,calincl->targetA(),calincl->targetZ());
+      proba = forceAbs(1,calincl->targetA(),calincl->targetZ(),calincl->bulletE(),bmaxt,proba_trans);
       
       standardRandom(&alea,&(hazard->igraine[4]));
@@ -7033 +7884 @@
       }
 
-      (*iarem) = int(calincl->f[0]) + bl3->ia1;
-      (*izrem) = int(calincl->f[1]) + int(iz1);
+      (*iarem) = int(calincl->targetA()) + bl3->ia1;
+      (*izrem) = int(calincl->targetZ()) + int(iz1);
       
-      del = std::sqrt(std::pow(((calincl->f[0]+1.)*fmpinc+calincl->f[2]),2)-pbeam2);
+      del = std::sqrt(std::pow(((calincl->targetA()+1.)*fmpinc+calincl->bulletE()),2)-pbeam2);
       
-      (*erecrem) = pbeam2/((calincl->f[0] + 1.0)*fmpinc + calincl->f[2] + del);
-
-      (*esrem) = calincl->f[2] + sep - (*erecrem);
+      (*erecrem) = pbeam2/((calincl->targetA() + 1.0)*fmpinc + calincl->bulletE() + del);
+
+      (*esrem) = calincl->bulletE() + sep - (*erecrem);
 
       (*alrem) = 0.00001;
@@ -7108 +7959 @@
     return dp0;
   }
-  G4double rp = radius(ap);                                                     
-  G4double rt = radius(at);                                                     
+  G4double rp = radius(G4int(ap));                                                     
+  G4double rt = radius(G4int(at));                                                     
   G4double vp = (dp1 + dpth)*dppi*std::pow(rp,3);                                          
   G4double vt = (dp1 + dpth)*dppi*std::pow(rt,3);                                          
@@ -7310 +8161 @@
 {
   // Gaussian random number generator
-  
+  G4double xg, xxg, xshuf;
+ gaussianRandom1:
+  xg=0.;
+  for(G4int i = 1; i <= 12; i++) { //do 6 j=1,12
+    standardRandom(&xxg,&(hazard->ial));
+    //  gaussianRandom6:
+    xg=xg+xxg;
+  }
+  xg=xg-6.;
+  if (xg*xg > 9.0) goto gaussianRandom1;
+  standardRandom(&xshuf,&(hazard->igraine[10]));
+  if (xshuf > 0.5) standardRandom(&xxg,&(hazard->ial));
+  // gaussianRandom2:
+  (*rndm) = xg;
+  return;
+    /*  
   G4double tempRandom = 0.0, random = 0.0, randomShuffle = 0.0;
  
@@ -7332 +8198 @@
 
   (*rndm) = random;
+    */
 }
 
@@ -7346 +8213 @@
 }
 
-G4double G4Incl::radius(G4double A)
+G4double G4Incl::radius(G4int A)
 {                                               
   const G4double dp1 = 1.0, dp3 = 3.0; 
@@ -7360 +8227 @@
              
   G4double fact = std::sqrt(dp5/dp3);                                                
-  G4int ia = int(std::floor(A+0.4));                                                        
+  //  G4int ia = int(std::floor(A+0.4));                                                        
+  G4int ia = A;
   G4double result = fact * (0.84 * std::pow(A,dpth) + 0.55);                               
   for(G4int i = 0; i < naSize; i++) {
@@ -7609 +8477 @@
 }
 
+// C****************************************************************************
+
+void G4Incl::projo_spec(G4int, G4int ips,
+                   G4double, G4double, G4double) // max_a_proj is a global #define'd as a global constant
+{
+  //      subroutine projo_spec(ia1,ips,max_a_proj,
+  //     sn_projspec,fmpinc,pinc,tlab)
+  //      common/bl1/p1(300),p2(300),p3(300),eps(300),ind1(300),ind2(300),ta
+  // spectators of the projectile treated here (nucleus for fermi-breakup):
+  // ips is the number of spect, their label is in n_projspec(ips).
+  // e,p are outside the target potential for these nucleons.
+  // la transmission de max_a_proj=60 ne marche pas....explicite ici!
+  //      integer a_projspec,z_projspec,n_projspec(60)
+  //      real*4 m_projspec
+  //      common/proj_spect/a_projspec,z_projspec,ex_projspec,t_projspec,
+  //     s p1_projspec,p2_projspec,p3_projspec,m_projspec
+  ///      common/quatvectprojo/eps_c(60),p3_c(60),
+  //     s  p1_s(60),p2_s(60),p3_s(60),
+  //     s  t_c(60)
+  // (eps_c,p1_s,p2_s,p3_s,eps_c used to store the kinematics of nucleons
+  // for composit projectiles before entering the potential)
+  //      common/fermi/tf,pf,pf2     
+
+  //  G4double ttrou[21];
+  G4double fja = 0.0;
+  G4double fjz = 0.0;
+  //  G4double ex2;
+  //  G4double rmmass;
+  //  const G4double um = 931.20793;
+
+// c verifs: ground state
+// c      es=0.
+// c      p1s=0.
+// c      p2s=0.
+// c      p3s=0.
+// c      gsmass=0.
+// c      do iloc=1,ia1
+// c      es=es+eps_c(iloc)
+// c      p1s=p1s+p1_s(iloc)
+// c      p2s=p2s+p2_s(iloc)
+// c      p3s=p3s+p3_s(iloc)
+// c      gsmass=gsmass+sqrt(eps_c(iloc)**2
+// c     s        -p1_s(iloc)**2-p2_s(iloc)**2-p3_s(iloc)**2)
+// c      enddo
+// c      sgs=sqrt(es**2-p1s**2-p2s**2-p3s**2)
+// c      write(6,*) 'sqs, gsmass :',sgs,gsmass
+// c      write(6,*) 'es,p1s,p2s,p3s:',es,p1s,p2s,p3s
+// c end verifs
+      
+//  G4cout <<"proj_spect called with: " << G4endl;
+//  G4cout <<"ips = " << ips << G4endl;
+  ps->a_projspec=0;
+  ps->z_projspec=0;
+  ps->t_projspec=0.0;
+  G4double e_spec=0.;
+    G4double p1_spec=0.;
+    G4double p2_spec=0.;
+    G4double p3_spec=0.;
+    G4double gs_mass=0.;
+
+  if(ips > 1) {
+    for(G4int i_spec=1; i_spec <= ips; i_spec++) {
+      G4int i_c = ps->n_projspec[i_spec];
+      e_spec = e_spec+bl1->eps[i_c];
+      p1_spec = p1_spec+bl1->p1[i_c];
+      p2_spec = p2_spec+bl1->p2[i_c];
+      p3_spec = p3_spec+bl1->p3[i_c];
+      G4double pi_spec2= std::pow(bl1->p1[i_c],2) + std::pow(bl1->p2[i_c],2) + std::pow(bl1->p3[i_c],2);
+      // current ground state mass:
+      gs_mass=gs_mass+std::sqrt(std::pow(bl1->eps[i_c],2)-pi_spec2);
+      if(bl1->ind2[i_c] == 1) ps->z_projspec=ps->z_projspec+1;
+    } // 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);
+
+    // no projectile spectator if a>=4 and a=z or a=n (no dresner breakup)
+    if(ips >= 4 && (ps->z_projspec == ips || ps->z_projspec == 0)) {
+      for(G4int i_spec = 1; i_spec <= ips; i_spec++) { //do i_spec=1,ips
+        G4int i_c = ps->n_projspec[i_spec];
+        // single nucleon forced on shell! momentum kept.
+        // single nucleon forced on shell! energy kept.
+        G4double arg=(std::pow(bl1->eps[i_c],2)-std::pow(938.2796,2))/
+          (std::pow(bl1->p1[i_c],2)+std::pow(bl1->p2[i_c], 2)+std::pow(bl1->p3[i_c],2));
+        G4double coef=0.0;
+        if(arg < 0.) { // Single nucleon is forced on shell. We have
+                       // to choose between energy and momentum
+                       // conservation... so we choose energy:
+          coef = 1.0;
+          bl1->eps[i_c]=std::sqrt(std::pow(bl1->p1[i_c],2) + std::pow(bl1->p2[i_c], 2) + 
+                                  std::pow(bl1->p3[i_c], 2) + std::pow(938.2796,2));
+        } else {
+          coef=std::sqrt(arg);
+          bl1->p1[i_c]=bl1->p1[i_c]*coef;
+          bl1->p2[i_c]=bl1->p2[i_c]*coef;
+          bl1->p3[i_c]=bl1->p3[i_c]*coef;
+        }
+      } // enddo
+      return;
+    } // endif
+        
+    // excitation energy from particle-hole excitation:
+    //   seconde methode de calcul ex4:
+    for(G4int i = 1; i <= bl3->ia1; i++) { // do i=1,ia1
+      G4int nb_p=i;
+      //      call ordered(t_c(i),nb_p,tpart)
+      ordered(qvp->t_c[i], nb_p);
+    } // enddo
+    G4double ex4=0.;
+    for(G4int i = 1; i <= ips; i++) { //do i=1,ips
+      G4int j=ps->n_projspec[i];
+      ex4=ex4+qvp->t_c[j]-ps->tab[i]; //ex4=ex4+qvp->t_c[j]-ps->tpart[i];
+    } // enddo
+
+    ps->ex_projspec =ex4;
+
+    // coherent rest mass:
+    ps->m_projspec = s_spec - ps->ex_projspec;
+
+    ps->a_projspec=ips;
+
+    if((ps->ex_projspec < 0.) || (ps->m_projspec < 0.)) {
+      G4cout <<"ex_projspec is negative,a,z: " << ps->ex_projspec << " " << fja << " " << fjz << G4endl;
+      G4cout <<"...or m is negative,a,z:" << ps->m_projspec << " " << fja << " " << fjz << G4endl;
+    } // endif
+
+    ps->t_projspec = e_spec-s_spec;
+    ps->p1_projspec = p1_spec;                  
+    ps->p2_projspec = p2_spec;
+    ps->p3_projspec = p3_spec;
+      
+// c check module...      
+// c      write(6,*) 'projectile spectator a,z:',a_projspec,z_projspec
+// c      write(6,*) 'm,t,e*:',m_projspec,t_projspec,ex_projspec
+// c      write(6,*) 'p1,p2,p3:',p1_projspec,p2_projspec,p3_projspec
+// c      write(6,*) '   constituants:',ips
+// c      do i_spec=1,ips
+// c            i_c = n_projspec(i_spec)
+// c      write(6,*) i_c,eps(i_c),p1(i_c),p2(i_c),p3(i_c)       
+// c      enddo
+// c ...end check module 
+
+    return;
+  } else if(ips == 1) { // then
+    ps->a_projspec = 0; // number of nucleons in the spectator nucleus
+    G4int i_c = ps->n_projspec[1];
+    //c single nucleon forced on shell! momentum kept.
+    //c      eps(i_c)=sqrt(p1(i_c)**2+p2(i_c)**2+p3(i_c)**2+938.2796**2)     
+    //c single nucleon forced on shell! energy kept.
+    G4double arg=(std::pow(bl1->eps[i_c],2)-std::pow(938.2796,2))/
+      (std::pow(bl1->p1[i_c],2)+std::pow(bl1->p2[i_c],2)+std::pow(bl1->p3[i_c],2));
+    G4double coef=0.0; // Initialize the coefficieant
+    if(arg < 0.) {
+      G4cout <<"off shell problems in projo_spec, forced coef = 1.0 to recover from the error" << G4endl;
+      coef = 1.0;
+    } else {
+      coef = std::sqrt(arg);
+    }
+
+    bl1->p1[i_c]=bl1->p1[i_c]*coef;
+    bl1->p2[i_c]=bl1->p2[i_c]*coef;
+    bl1->p3[i_c]=bl1->p3[i_c]*coef;
+  } //endif
+
+  ps->a_projspec = 0; // number of nucleons in the spectator nucleus
+  return;
+}
+      
+void G4Incl::ordered(G4double t, G4int nb)
+{
+  //      SUBROUTINE ORDERED(T,nb,TAB)
+  //C This subroutine will ordered in TAB() from min TAB(1) to max
+  //C TAB(nb), nb values of T sent successively. 
+  //C CALL ORDERED (T1,1), CALL ORDERED (T2,2),..CALL ORDERED (Tnb,nb) 
+  //      DIMENSION TAB(1)
+
+  //  G4cout <<"Ordered was called" << G4endl;
+  if(nb == 1) {
+    ps->tab[1]=t;
+    return;
+  } // endif
+
+  G4int i=0;
+
+ ordered100:
+  i=i+1;
+  if(t < ps->tab[i]) {
+    //    for(G4int l = nb-1; l != i; l--) { //do l=nb-1,i,-1
+    for(G4int l = nb-1; l >= i; l--) { //do l=nb-1,i,-1
+      //      G4cout <<"Ordered: L = " << l << " nb = " << nb <<  " i = " << i << G4endl;
+      ps->tab[l+1]=ps->tab[l];
+    } // enddo
+    ps->tab[i]=t;
+    return;
+  } else {
+    if(i <= nb-1) goto ordered100;
+    ps->tab[i]=t;
+    return;
+  } // endif
+  return;
+}
+
 // Utilities
 
@@ -7903 +8972 @@
   G4cout <<"(particle " << index << G4endl;
   if(bl1->ind1[index] == 0) { // Nucleon
-    if(bl1->ind1[index] == 1) G4cout << "(quote proton)" << G4endl;
-    else if(bl1->ind1[index] == -1) G4cout << "(quote neutron)" << G4endl;
-    else G4cout << "(quote unidentified-nucleon)" << G4endl;
+    if(bl1->ind2[index] == 1) {
+      G4cout << "(quote proton)" << G4endl;
+    } else if(bl1->ind2[index] == -1) {
+      G4cout << "(quote neutron)" << G4endl;
+    } else {
+      G4cout << "(quote unidentified-nucleon-ind1-" << bl1->ind1[index] <<"-ind2-" << bl1->ind2[index] << G4endl;
+    }
   } else if(bl1->ind1[index] == 1) { // Delta
-    if(bl1->ind2[index] >= -2 && bl1->ind2[index] <= 2) G4cout << "(quote delta-"
-                                                               << bl1->ind2[index]
-                                                               << ")" << G4endl;
-    else G4cout << "(quote unidentified-delta)" << G4endl;
+    if(bl1->ind2[index] >= -2 && bl1->ind2[index] <= 2) {
+      G4cout << "(quote delta-"
+             << bl1->ind2[index]
+             << ")" << G4endl;
+    } else {
+      G4cout << "(quote unidentified-delta)" << G4endl;
+    }
   } else {
     G4cout <<"(quote unidentified-particle)" << G4endl;
@@ -7938 +9014 @@
   G4cout <<"(list ;; Particles at the beginning of the avatar" << G4endl;
   G4int ia = bl3->ia1 + bl3->ia2;
-  for(G4int i = 1; i != ia; ++i) { // do i=1,IA
+  for(G4int i = 1; i <= ia; ++i) { // do i=1,IA
     print_one_particle(i);
   }

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4InclAblaCascadeInterface.cc,v 1.13 2009/12/04 13:16:57 kaitanie Exp $ 
+// $Id: G4InclAblaCascadeInterface.cc,v 1.16 2010/10/29 06:48:43 gunter Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -34 +34 @@
 
 #include "G4InclAblaCascadeInterface.hh"
+#include "G4FermiBreakUp.hh"
 #include "math.h"
 #include "G4GenericIon.hh"
@@ -46 +47 @@
 
   varntp = new G4VarNtp();
-  calincl = new G4Calincl();
+  calincl = 0;
   ws = new G4Ws();
   mat = new G4Mat();
   incl = new G4Incl(hazard, calincl, ws, mat, varntp);
+  if(!getenv("G4INCLABLANOFERMIBREAKUP")) { // Use Fermi Break-up by default if it is NOT explicitly disabled
+    incl->setUseFermiBreakUp(true);
+  }
 
   verboseLevel = 0;
@@ -58 +62 @@
   delete hazard;
   delete varntp;
-  delete calincl;
   delete ws;
   delete mat;
@@ -68 +71 @@
   G4int maxTries = 200;
 
-  G4int particleI, n = 0;
-
+  G4int particleI;
   G4int bulletType = 0;
 
@@ -86 +88 @@
   }
 
-  // INCL4 needs the energy in units MeV
-  G4double bulletE = aTrack.GetKineticEnergy() * MeV;
-
 #ifdef DEBUGINCL
   G4cout <<"Bullet energy = " << bulletE / MeV << G4endl;
 #endif
-
-  G4double targetA = theNucleus.GetN();
-  G4double targetZ = theNucleus.GetZ();
 
   G4double eKin;
@@ -100 +96 @@
   G4DynamicParticle *cascadeParticle = 0;
   G4ParticleDefinition *aParticleDefinition = 0;
+  
+  G4FermiBreakUp *fermiBreakUp = new G4FermiBreakUp();
+  G4FragmentVector *theFermiBreakupResult = 0;
+  G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
 
   // INCL assumes the projectile particle is going in the direction of
@@ -113 +113 @@
   theResult.Clear(); // Make sure the output data structure is clean.
 
-  // 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() ==    G4Proton::Proton()    ) bulletType = proton;
-  if (aTrack.GetDefinition() ==   G4Neutron::Neutron()   ) bulletType = neutron;
-  if (aTrack.GetDefinition() ==  G4PionPlus::PionPlus()  ) bulletType = pionPlus;
-  if (aTrack.GetDefinition() == G4PionMinus::PionMinus() ) bulletType = pionMinus;
-  if (aTrack.GetDefinition() ==  G4PionZero::PionZero()  ) bulletType = pionZero;
+  calincl = new G4InclInput(aTrack, theNucleus, false);
+  incl->setInput(calincl);
+
+  G4InclInput::printProjectileTargetInfo(aTrack, theNucleus);
+  calincl->printInfo();
 
 #ifdef DEBUGINCL
@@ -139 +134 @@
 #endif
 
- for(int i = 0; i < 15; i++) {
-   calincl->f[i] = 0.0; // Initialize INCL input data
- }
-
   // Check wheter the input is acceptable.
-  if((bulletType != 0) && ((targetA != 1) && (targetZ != 1))) {
-    calincl->f[0] = targetA;    // Target mass number
-    calincl->f[1] = targetZ;    // Charge number
-    calincl->f[6] = bulletType; // Type
-    calincl->f[2] = bulletE;    // Energy [MeV]
-    calincl->f[5] = 1.0;        // Time scaling
-    calincl->f[4] = 45.0;       // Nuclear potential
-
+  if((calincl->bulletType() != 0) && ((calincl->targetA() != 1) && (calincl->targetZ() != 1))) {
     ws->nosurf = -2;  // Nucleus surface, -2 = Woods-Saxon 
     ws->xfoisa = 8;
@@ -160 +144 @@
 
     mat->nbmat = 1;
-    mat->amat[0] = int(calincl->f[0]);
-    mat->zmat[0] = int(calincl->f[1]);
+    mat->amat[0] = int(calincl->targetA());
+    mat->zmat[0] = int(calincl->targetZ());
 
     incl->initIncl(true);
@@ -170 +154 @@
         G4cout <<"G4InclAblaCascadeInterface: Try number = " << nTries << G4endl; 
       }
-      incl->processEventInclAbla(eventNumber);
+      incl->processEventInclAbla(calincl, eventNumber);
 
       if(verboseLevel > 1) {
@@ -181 +165 @@
        * Diagnostic output
        */
-      G4cout <<"G4InclAblaCascadeInterface: Bullet type: " << bulletType << G4endl;
-      G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << bulletE << " MeV" << G4endl;
-
-      G4cout <<"G4InclAblaCascadeInterface: Target A:  " << targetA << G4endl;
-      G4cout <<"G4InclAblaCascadeInterface: Target Z:  " << targetZ << G4endl;
+      G4cout <<"G4InclAblaCascadeInterface: Bullet type: " << calincl->bulletType() << G4endl;
+      G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl;
+
+      G4cout <<"G4InclAblaCascadeInterface: Target A:  " << calincl->targetA() << G4endl;
+      G4cout <<"G4InclAblaCascadeInterface: Target Z:  " << calincl->targetZ() << G4endl;
 
       if(verboseLevel > 3) {
-        diagdata <<"G4InclAblaCascadeInterface: Bullet type: " << bulletType << G4endl;
-        diagdata <<"G4InclAblaCascadeInterface: Bullet energy: " << bulletE << " MeV" << G4endl;
+        diagdata <<"G4InclAblaCascadeInterface: Bullet type: " << calincl->bulletType() << G4endl;
+        diagdata <<"G4InclAblaCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl;
         
-        diagdata <<"G4InclAblaCascadeInterface: Target A:  " << targetA << G4endl;
-        diagdata <<"G4InclAblaCascadeInterface: Target Z:  " << targetZ << G4endl;
-      }
-
-      for(particleI = 0; particleI < varntp->ntrack; particleI++) {
-        G4cout << n                       << " " << calincl->f[6]             << " " << calincl->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                       << " " << calincl->f[6]             << " " << calincl->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 <<"G4InclAblaCascadeInterface: Target A:  " << calincl->targetA() << G4endl;
+        diagdata <<"G4InclAblaCascadeInterface: Target Z:  " << calincl->targetZ() << G4endl;
+      }
+
     }
 
@@ -227 +191 @@
       theResult.SetStatusChange(stopAndKill);
       
-      if(bulletType == proton) {
-        aParticleDefinition = G4Proton::ProtonDefinition();
-      } else if(bulletType == neutron) {
-        aParticleDefinition = G4Neutron::NeutronDefinition();
-      } else if(bulletType == pionPlus) {
-        aParticleDefinition = G4PionPlus::PionPlusDefinition();
-      } else if(bulletType == pionZero) {
-        aParticleDefinition = G4PionZero::PionZeroDefinition();
-      } else if(bulletType == pionMinus) {
-        aParticleDefinition = G4PionMinus::PionMinusDefinition();
-      } else { // Projectile was not regognized
-        aParticleDefinition = 0;
-      }
+      G4int bulletType = calincl->bulletType();
+      aParticleDefinition = G4InclInput::getParticleDefinition(bulletType);
 
       if(aParticleDefinition != 0) {
@@ -328 +281 @@
       if((varntp->avv[particleI] > 1) && (varntp->zvv[particleI] >= 1)) { // Nucleus fragment
         G4ParticleDefinition * aIonDef = 0;
-        G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
 
         G4int A = G4int(varntp->avv[particleI]);
@@ -404 +356 @@
         }
       }
+    }
+
+    // Finally do Fermi break-up if needed
+    if(varntp->needsFermiBreakup) {
+      //      baryonNumberBalanceInINCL -= varntp->massini;
+      //      chargeNumberBalanceInINCL -= 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);
+      theFermiBreakupResult = fermiBreakUp->BreakItUp(theCascadeRemnant);
+      if(theFermiBreakupResult != 0) {
+      G4FragmentVector::iterator fragment;
+      for(fragment = theFermiBreakupResult->begin(); fragment != theFermiBreakupResult->end(); fragment++) {
+        G4ParticleDefinition *theFragmentDefinition = 0;
+        if((*fragment)->GetA_asInt() == 1 && (*fragment)->GetZ_asInt() == 0) { // Neutron
+          theFragmentDefinition = G4Neutron::NeutronDefinition();
+        } else if ((*fragment)->GetA_asInt() == 1 && (*fragment)->GetZ_asInt() == 1) {
+          theFragmentDefinition = G4Proton::ProtonDefinition();
+        } else {
+          theFragmentDefinition = theTableOfParticles->GetIon((*fragment)->GetZ_asInt(), (*fragment)->GetA_asInt(), (*fragment)->GetExcitationEnergy());
+        }
+
+        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 <<"G4InclAblaCascadeInterface: Error. Fragment produced by Fermi break-up does not exist." << G4endl;
+          G4cout <<"Resulting fragment: " << G4endl;
+          G4cout <<" Z = " << (*fragment)->GetZ_asInt() << G4endl;
+          G4cout <<" A = " << (*fragment)->GetA_asInt() << G4endl;
+          G4cout <<" Excitation : " << (*fragment)->GetExcitationEnergy() / MeV << " MeV" << G4endl;
+          G4cout <<" momentum = " << (*fragment)->GetMomentum().mag() / MeV << " MeV" << G4endl;
+        }
+      }
+      if(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) {
+        G4cout <<"Baryon number balance after remnant nucleus Fermi break-up: " << baryonNumberBalance << G4endl;
+      }
+      if(chargeBalance != 0) {
+        G4cout <<"Charge balance after remnant nucleus Fermi break-up: " << chargeBalance << G4endl;
+      }
+    }
     }
 
@@ -449 +500 @@
     theResult.SetStatusChange(stopAndKill);
 
-    G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();
     cascadeParticle = new G4DynamicParticle(theTableOfParticles->FindParticle(aTrack.GetDefinition()), aTrack.Get4Momentum());
 
@@ -472 +522 @@
     }
 
-    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;
       }
     }
@@ -500 +550 @@
   }
 
+  delete fermiBreakUp;
+  delete calincl;
+  calincl = 0;
   return &theResult;
 } 

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4InclAblaLightIonInterface.cc,v 1.11 2009/12/04 13:16:57 kaitanie Exp $ 
+// $Id: G4InclAblaLightIonInterface.cc,v 1.14 2010/10/29 06:48:43 gunter Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -31 +31 @@
 // Aatos Heikkinen, HIP (project coordination)
 
+#include <vector>
+
 #include "G4InclAblaLightIonInterface.hh"
+#include "G4FermiBreakUp.hh"
 #include "math.h"
 #include "G4GenericIon.hh"
@@ -44 +47 @@
 
   varntp = new G4VarNtp();
-  calincl = new G4Calincl();
+  calincl = 0;
   ws = new G4Ws();
   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;
+  }
 }
 
@@ -64 +76 @@
 G4HadFinalState* G4InclAblaLightIonInterface::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
 {
+  //  const G4bool useFermiBreakup = false;
   G4int maxTries = 200;
 
-  G4int particleI, n = 0;
-
-  G4int bulletType = 0;
-
-  // Print diagnostic messages: 0 = silent, 1 and 2 = verbose
-  verboseLevel = 0;
+  G4int particleI;
+
+  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 << " >>> G4InclAblaLightIonInterface::ApplyYourself called" << G4endl;
@@ -84 +103 @@
   }
 
-  // INCL4 needs the energy in units MeV
-  G4double bulletE = aTrack.GetKineticEnergy() / MeV;
-
-  G4double targetA = theNucleus.GetN();
-  G4double targetZ = theNucleus.GetZ();
+  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)) {
+    G4ParticleDefinition *oldTargetDef = theTableOfParticles->GetIon(theNucleus.GetA_asInt(), theNucleus.GetZ_asInt(), 0.0);
+    const G4ParticleDefinition *oldProjectileDef = aTrack.GetDefinition();
+
+    G4int oldTargetA = oldTargetDef->GetAtomicMass();
+    G4int newTargetA = oldProjectileDef->GetAtomicMass();
+    G4int newTargetZ = oldProjectileDef->GetAtomicNumber();
+
+    if(newTargetA > 0 && newTargetZ > 0 && oldTargetDef != 0 && oldProjectileDef != 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;
@@ -105 +157 @@
   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;
+  */
+
+  G4FermiBreakUp *fermiBreakUp = new G4FermiBreakUp();
+  G4FragmentVector *theSpectatorFermiBreakupResult = 0;
+  G4FragmentVector *theFermiBreakupResult = 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;
-
-  for(int i = 0; i < 15; i++) {
-    calincl->f[i] = 0.0; // Initialize INCL input data
+                           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))) {
-    calincl->f[0] = targetA;    // Target mass number
-    calincl->f[1] = targetZ;    // Charge number
-    calincl->f[6] = bulletType; // Type
-    calincl->f[2] = bulletE;    // Energy [MeV]
-    calincl->f[5] = 1.0;        // Time scaling
-    calincl->f[4] = 45.0;       // Nuclear potential
-
+  if((bulletType != 0) && ((calincl->targetA() != 1) && (calincl->targetZ() != 1))) {
     ws->nosurf = -2;  // Nucleus surface, -2 = Woods-Saxon 
     ws->xfoisa = 8;
@@ -138 +224 @@
 
     mat->nbmat = 1;
-    mat->amat[0] = int(calincl->f[0]);
-    mat->zmat[0] = int(calincl->f[1]);
-
+    mat->amat[0] = int(calincl->targetA());
+    mat->zmat[0] = int(calincl->targetA());
+
+    incl->setInput(calincl);
     incl->initIncl(true);
 
@@ -148 +235 @@
         G4cout <<"G4InclAblaLightIonInterface: Try number = " << nTries << G4endl; 
       }
-      incl->processEventInclAbla(eventNumber);
+      incl->processEventInclAbla(calincl, eventNumber);
 
       if(verboseLevel > 1) {
@@ -159 +246 @@
        * Diagnostic output
        */
-      G4cout <<"G4InclAblaLightIonInterface: Bullet type: " << bulletType << G4endl;
-      G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << bulletE << " MeV" << G4endl;
-
-      G4cout <<"G4InclAblaLightIonInterface: Target A:  " << targetA << G4endl;
-      G4cout <<"G4InclAblaLightIonInterface: Target Z:  " << targetZ << G4endl;
+      G4cout <<"G4InclAblaLightIonInterface: 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 <<"G4InclAblaLightIonInterface: Target A:  " << calincl->targetA() << G4endl;
+      G4cout <<"G4InclAblaLightIonInterface: Target Z:  " << calincl->targetZ() << G4endl;
 
       if(verboseLevel > 3) {
-        diagdata <<"G4InclAblaLightIonInterface: Bullet type: " << bulletType << G4endl;
-        diagdata <<"G4InclAblaLightIonInterface: Bullet energy: " << bulletE << " MeV" << G4endl;
+        diagdata <<"G4InclAblaLightIonInterface: Bullet type: " << calincl->bulletType() << G4endl;
+        diagdata <<"G4InclAblaLightIonInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl;
         
-        diagdata <<"G4InclAblaLightIonInterface: Target A:  " << targetA << G4endl;
-        diagdata <<"G4InclAblaLightIonInterface: Target Z:  " << targetZ << G4endl;
-      }
-
-      for(particleI = 0; particleI < varntp->ntrack; particleI++) {
-        G4cout << n                       << " " << calincl->f[6]             << " " << calincl->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                       << " " << calincl->f[6]             << " " << calincl->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 <<"G4InclAblaLightIonInterface: Target A:  " << calincl->targetA() << G4endl;
+        diagdata <<"G4InclAblaLightIonInterface: Target Z:  " << calincl->targetZ() << G4endl;
       }
     }
@@ -231 +301 @@
         cascadeParticle->SetDefinition(aParticleDefinition);
         cascadeParticle->Set4Momentum(aTrack.Get4Momentum());
-        theResult.AddSecondary(cascadeParticle);
+        result.push_back(cascadeParticle);
       }
     }
@@ -239 +309 @@
     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;
@@ -262 +334 @@
           new G4DynamicParticle(G4Proton::ProtonDefinition(), momDirection, eKin);
         particleIdentified++;
+        baryonNumberBalanceInINCL -= 1;
+        chargeNumberBalanceInINCL -= 1;
       }
 
@@ -268 +342 @@
           new G4DynamicParticle(G4Neutron::NeutronDefinition(), momDirection, eKin);
         particleIdentified++;
+        baryonNumberBalanceInINCL -= 1;
       }
 
@@ -274 +349 @@
           new G4DynamicParticle(G4PionPlus::PionPlusDefinition(), momDirection, eKin);
         particleIdentified++;
+        chargeNumberBalanceInINCL -= 1;
       }
 
@@ -280 +356 @@
           new G4DynamicParticle(G4PionZero::PionZeroDefinition(), momDirection, eKin);
         particleIdentified++;
+        chargeNumberBalanceInINCL -= 0;
       }
 
@@ -286 +363 @@
           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]);
@@ -313 +390 @@
             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.
@@ -332 +411 @@
     }
 
+    // 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);
+      theSpectatorFermiBreakupResult = fermiBreakUp->BreakItUp(theSpectatorNucleus);
+      if(theSpectatorFermiBreakupResult != 0) {
+      G4FragmentVector::iterator fragment;
+      for(fragment = theSpectatorFermiBreakupResult->begin(); fragment != theSpectatorFermiBreakupResult->end(); fragment++) {
+        G4ParticleDefinition *theFragmentDefinition = 0;
+        if((*fragment)->GetA_asInt() == 1 && (*fragment)->GetZ_asInt() == 0) { // Neutron
+          theFragmentDefinition = G4Neutron::NeutronDefinition();
+        } else if ((*fragment)->GetA_asInt() == 1 && (*fragment)->GetZ_asInt() == 1) {
+          theFragmentDefinition = G4Proton::ProtonDefinition();
+        } else {
+          theFragmentDefinition = theTableOfParticles->GetIon((*fragment)->GetZ_asInt(), (*fragment)->GetA_asInt(), (*fragment)->GetExcitationEnergy());
+        }
+        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;
+          }
+          result.push_back(theFragment);
+        } else {
+          G4cout <<"G4InclAblaCascadeInterface: Error. Fragment produced by Fermi break-up does not exist." 
+                 << G4endl;
+          G4cout <<"Resulting fragment: " << G4endl;
+          G4cout <<" Z = " << (*fragment)->GetZ_asInt() << G4endl;
+          G4cout <<" A = " << (*fragment)->GetA_asInt() << G4endl;
+          G4cout <<" Excitation : " << (*fragment)->GetExcitationEnergy() / MeV << " MeV" << G4endl;
+          G4cout <<" momentum = " << (*fragment)->GetMomentum().mag() / MeV << " MeV" << G4endl;
+        }
+      }
+      if(std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 
+        G4cout <<"Four-momentum balance after spectator 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) {
+        G4cout <<"Event " << eventNumber << ": Baryon number balance after spectator nucleus Fermi break-up: " << baryonNumberBalance << G4endl;
+      }
+      if(chargeBalance != 0) {
+        G4cout <<"Event " << eventNumber <<": Charge balance after spectator nucleus Fermi break-up: " << chargeBalance << G4endl;
+      }
+    }
+  }
+
+    // Finally do Fermi break-up if needed
+    if(varntp->needsFermiBreakup && 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);
+      theFermiBreakupResult = fermiBreakUp->BreakItUp(theCascadeRemnant);
+      if(theFermiBreakupResult != 0) {
+      G4FragmentVector::iterator fragment;
+      for(fragment = theFermiBreakupResult->begin(); fragment != theFermiBreakupResult->end(); fragment++) {
+        G4ParticleDefinition *theFragmentDefinition = 0;
+        if((*fragment)->GetA_asInt() == 1 && (*fragment)->GetZ_asInt() == 0) { // Neutron
+          theFragmentDefinition = G4Neutron::NeutronDefinition();
+        } else if ((*fragment)->GetA_asInt() == 1 && (*fragment)->GetZ_asInt() == 1) {
+          theFragmentDefinition = G4Proton::ProtonDefinition();
+        } else {
+          theFragmentDefinition = theTableOfParticles->GetIon((*fragment)->GetZ_asInt(), (*fragment)->GetA_asInt(), (*fragment)->GetExcitationEnergy());
+        }
+
+        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;
+          }
+          result.push_back(theFragment);
+        } else {
+          G4cout <<"G4InclAblaCascadeInterface: Error. Fragment produced by Fermi break-up does not exist." << G4endl;
+          G4cout <<"Resulting fragment: " << G4endl;
+          G4cout <<" Z = " << (*fragment)->GetZ_asInt() << G4endl;
+          G4cout <<" A = " << (*fragment)->GetA_asInt() << G4endl;
+          G4cout <<" Excitation : " << (*fragment)->GetExcitationEnergy() / MeV << " MeV" << G4endl;
+          G4cout <<" momentum = " << (*fragment)->GetMomentum().mag() / MeV << " MeV" << G4endl;
+        }
+      }
+      if(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) {
+        G4cout <<"Baryon number balance after remnant nucleus Fermi break-up: " << baryonNumberBalance << G4endl;
+      }
+      if(chargeBalance != 0) {
+        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 <<": G4InclAblaLightIonInterface: 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 <<": G4InclAblaLightIonInterface: Charge number conservation problem in INCL detected!" << G4endl;
+      G4cout <<"Event " << eventNumber <<": Charge number balance: " << chargeNumberBalanceInINCL << G4endl;
+    }
   }
   /**
@@ -344 +633 @@
     cascadeParticle = new G4DynamicParticle(theTableOfParticles->FindParticle(aTrack.GetDefinition()), aTrack.Get4Momentum());
 
-    theResult.AddSecondary(cascadeParticle);
+    result.push_back(cascadeParticle);
 
     if(verboseLevel > 1) {
@@ -364 +653 @@
     }
 
-    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;
       }
     }
@@ -392 +681 @@
   }
 
+  // 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;
 } 

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4InclCascadeInterface.cc,v 1.10 2007/12/10 16:32:02 gunter Exp $ 
+// $Id: G4InclCascadeInterface.cc,v 1.12 2010/10/26 02:47:59 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -45 +45 @@
 
   varntp = new G4VarNtp();
-  calincl = new G4Calincl();
+  inclInput = 0;
   ws = new G4Ws();
   mat = new G4Mat();
-  incl = new G4Incl(hazard, calincl, ws, mat, varntp);
+  incl = new G4Incl(hazard, inclInput, ws, mat, varntp);
 
   verboseLevel = 0;
@@ -57 +57 @@
   delete hazard;
   delete varntp;
-  delete calincl;
+  delete inclInput;
   delete ws;
   delete mat;
@@ -67 +67 @@
   G4int maxTries = 200;
 
-  G4int particleI, n = 0;
-
-  G4int bulletType = 0;
+  G4int particleI;
 
   // Print diagnostic messages: 0 = silent, 1 and 2 = verbose
@@ -85 +83 @@
   }
 
-  // INCL4 needs the energy in units MeV
-  G4double bulletE = aTrack.GetKineticEnergy() * MeV;
-
 #ifdef DEBUGINCL
   G4cout <<"Bullet energy = " << bulletE / MeV << G4endl;
 #endif
-
-  G4double targetA = theNucleus.GetN();
-  G4double targetZ = theNucleus.GetZ();
 
   G4double eKin;
@@ -110 +102 @@
   G4LorentzRotation toLabFrame = toZ.inverse();
 
+  inclInput = new G4InclInput(aTrack, theNucleus, false);
+
   theResult.Clear(); // Make sure the output data structure is clean.
-
-  // 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() ==    G4Proton::Proton()    ) bulletType = proton;
-  if (aTrack.GetDefinition() ==   G4Neutron::Neutron()   ) bulletType = neutron;
-  if (aTrack.GetDefinition() ==  G4PionPlus::PionPlus()  ) bulletType = pionPlus;
-  if (aTrack.GetDefinition() == G4PionMinus::PionMinus() ) bulletType = pionMinus;
-  if (aTrack.GetDefinition() ==  G4PionZero::PionZero()  ) bulletType = pionZero;
 
 #ifdef DEBUGINCL
@@ -137 +120 @@
 #endif
 
-  for(int i = 0; i < 15; i++) {
-    calincl->f[i] = 0.0; // Initialize INCL input data
-  }
-
   // Check wheter the input is acceptable.
-  if((bulletType != 0) && ((targetA != 1) && (targetZ != 1))) {
-    calincl->f[0] = targetA;    // Target mass number
-    calincl->f[1] = targetZ;    // Charge number
-    calincl->f[6] = bulletType; // Type
-    calincl->f[2] = bulletE;    // Energy [MeV]
-    calincl->f[5] = 1.0;        // Time scaling
-    calincl->f[4] = 45.0;       // Nuclear potential
-
+  if((inclInput->bulletType() != 0) && ((inclInput->targetA() != 1) && (inclInput->targetZ() != 1))) {
     ws->nosurf = -2;  // Nucleus surface, -2 = Woods-Saxon 
     ws->xfoisa = 8;
@@ -158 +130 @@
 
     mat->nbmat = 1;
-    mat->amat[0] = int(calincl->f[0]);
-    mat->zmat[0] = int(calincl->f[1]);
+    mat->amat[0] = int(inclInput->targetA());
+    mat->zmat[0] = int(inclInput->targetZ());
 
     incl->initIncl(true);
@@ -168 +140 @@
         G4cout <<"G4InclCascadeInterface: Try number = " << nTries << G4endl; 
       }
-      incl->processEventIncl();
+      incl->processEventIncl(inclInput);
 
       if(verboseLevel > 1) {
@@ -179 +151 @@
        * Diagnostic output
        */
-      G4cout <<"G4InclCascadeInterface: Bullet type: " << bulletType << G4endl;
-      G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << bulletE << " MeV" << G4endl;
-
-      G4cout <<"G4InclCascadeInterface: Target A:  " << targetA << G4endl;
-      G4cout <<"G4InclCascadeInterface: Target Z:  " << targetZ << G4endl;
+      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;
 
       if(verboseLevel > 3) {
-        diagdata <<"G4InclCascadeInterface: Bullet type: " << bulletType << G4endl;
-        diagdata <<"G4InclCascadeInterface: Bullet energy: " << bulletE << " MeV" << G4endl;
+        diagdata <<"G4InclCascadeInterface: Bullet type: " << inclInput->bulletType() << G4endl;
+        diagdata <<"G4InclCascadeInterface: Bullet energy: " << inclInput->bulletE() << " MeV" << G4endl;
         
-        diagdata <<"G4InclCascadeInterface: Target A:  " << targetA << G4endl;
-        diagdata <<"G4InclCascadeInterface: Target Z:  " << targetZ << G4endl;
-      }
-
-      for(particleI = 0; particleI < varntp->ntrack; particleI++) {
-        G4cout << n                       << " " << calincl->f[6]             << " " << calincl->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                       << " " << calincl->f[6]             << " " << calincl->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:  " << 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;
+      //   }
+      // }
     }
 
@@ -224 +196 @@
 
       theResult.SetStatusChange(stopAndKill);
-      
-      if(bulletType == proton) {
-        aParticleDefinition = G4Proton::ProtonDefinition();
-      }
-      if(bulletType == neutron) {
-        aParticleDefinition = G4Neutron::NeutronDefinition();
-      }
-      if(bulletType == pionPlus) {
-        aParticleDefinition = G4PionPlus::PionPlusDefinition();
-      }
-      if(bulletType == pionZero) {
-        aParticleDefinition = G4PionZero::PionZeroDefinition();
-      }
-      if(bulletType == pionMinus) {
-        aParticleDefinition = G4PionMinus::PionMinusDefinition();
-      }
+      aParticleDefinition = G4InclInput::getParticleDefinition(inclInput->bulletType());
 
       cascadeParticle = new G4DynamicParticle();
@@ -434 +391 @@
     }
 
-    if(bulletType == 0) {
+    if(inclInput->bulletType() == 0) {
       if(verboseLevel > 1) {
         G4cout <<"G4InclCascadeInterface: Unknown bullet type" << G4endl;      
@@ -445 +402 @@
     }
 
-    if((targetA == 1) && (targetZ == 1)) { // Unsupported target
+    if((inclInput->targetA() == 1) && (inclInput->targetZ() == 1)) { // Unsupported target
       if(verboseLevel > 1) {
         G4cout <<"Unsupported target: " << G4endl;
-        G4cout <<"Target A: " << targetA << G4endl;
-        G4cout <<"TargetZ: " << targetZ << G4endl;
+        G4cout <<"Target A: " << inclInput->targetA() << G4endl;
+        G4cout <<"TargetZ: " << inclInput->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.
-      if(verboseLevel > 1) {
-        G4cout <<"Unsupported bullet energy: " << bulletE << " MeV. (Lower limit is 100 MeV)." << G4endl;
+        diagdata <<"Target A: " << inclInput->targetA() << G4endl;
+        diagdata <<"TargetZ: " << inclInput->targetZ() << G4endl;
+      }
+    }
+
+    if(inclInput->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 <<"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: " << inclInput->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;
       }
     }

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4InclLightIonInterface.cc,v 1.10 2007/12/10 16:32:07 gunter Exp $ 
+// $Id: G4InclLightIonInterface.cc,v 1.12 2010/10/26 02:47:59 kaitanie Exp $ 
 // Translation of INCL4.2/ABLA V3 
 // Pekka Kaitaniemi, HIP (translation)
@@ -43 +43 @@
 
   varntp = new G4VarNtp();
-  calincl = new G4Calincl();
+  calincl = 0;
   ws = new G4Ws();
   mat = new G4Mat();
@@ -65 +65 @@
   G4int maxTries = 200;
 
-  G4int particleI, n = 0;
-
   G4int bulletType = 0;
+  G4int particleI;
 
   // Print diagnostic messages: 0 = silent, 1 and 2 = verbose
@@ -86 +85 @@
   G4double bulletE = aTrack.GetKineticEnergy() / MeV;
 
-  G4double targetA = theNucleus.GetN();
-  G4double targetZ = theNucleus.GetZ();
+  G4int targetA = theNucleus.GetA_asInt();
+  G4int targetZ = theNucleus.GetZ_asInt();
 
   G4double eKin;
@@ -116 +115 @@
   if (aTrack.GetDefinition() == G4Alpha::Alpha()         ) bulletType = he4;
 
-  for(int i = 0; i < 15; i++) {
-    calincl->f[i] = 0.0; // Initialize INCL input data
-  }
+  calincl = new G4InclInput();
 
   // Check wheter the input is acceptable.
   if((bulletType != 0) && ((targetA != 1) && (targetZ != 1))) {
-    calincl->f[0] = targetA;    // Target mass number
-    calincl->f[1] = targetZ;    // Charge number
-    calincl->f[6] = bulletType; // Type
-    calincl->f[2] = bulletE;    // Energy [MeV]
-    calincl->f[5] = 1.0;        // Time scaling
-    calincl->f[4] = 45.0;       // Nuclear potential
-
     ws->nosurf = -2;  // Nucleus surface, -2 = Woods-Saxon 
     ws->xfoisa = 8;
@@ -137 +127 @@
 
     mat->nbmat = 1;
-    mat->amat[0] = int(calincl->f[0]);
-    mat->zmat[0] = int(calincl->f[1]);
+    mat->amat[0] = int(calincl->targetA());
+    mat->zmat[0] = int(calincl->targetZ());
 
     incl->initIncl(true);
@@ -147 +137 @@
         G4cout <<"G4InclLightIonInterface: Try number = " << nTries << G4endl; 
       }
-      incl->processEventIncl();
+      incl->processEventIncl(calincl);
 
       if(verboseLevel > 1) {
@@ -170 +160 @@
         diagdata <<"G4InclLightIonInterface: Target A:  " << targetA << G4endl;
         diagdata <<"G4InclLightIonInterface: Target Z:  " << targetZ << G4endl;
-      }
-
-      for(particleI = 0; particleI < varntp->ntrack; particleI++) {
-        G4cout << n                       << " " << calincl->f[6]             << " " << calincl->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                       << " " << calincl->f[6]             << " " << calincl->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;
-        }
       }
     }