Index: trunk/source/processes/hadronic/models/theo_high_energy/History
===================================================================
--- trunk/source/processes/hadronic/models/theo_high_energy/History	(revision 1007)
+++ trunk/source/processes/hadronic/models/theo_high_energy/History	(revision 1055)
@@ -14,4 +14,8 @@
      * Please list in reverse chronological order (last date on top)
      ---------------------------------------------------------------
+
+29-Apr-2009 M.Kosov         had-theo-HE-V09-02-02
+- Quasi Elastic part is made more used in QGSC at low energies to improve
+  the performance of QGSC at low energies
 
 8-Apr-2008  G.Folger        had-theo-HE-V09-01-01
Index: trunk/source/processes/hadronic/models/theo_high_energy/src/G4QuasiElasticChannel.cc
===================================================================
--- trunk/source/processes/hadronic/models/theo_high_energy/src/G4QuasiElasticChannel.cc	(revision 1007)
+++ trunk/source/processes/hadronic/models/theo_high_energy/src/G4QuasiElasticChannel.cc	(revision 1055)
@@ -25,9 +25,10 @@
 //
 //
-// $Id: G4QuasiElasticChannel.cc,v 1.4 2008/04/24 13:26:19 gunter Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4QuasiElasticChannel.cc,v 1.7 2009/04/09 08:28:42 mkossov Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 
 // Author : Gunter Folger March 2007
+// Modified by Mikhail Kossov. Apr2009, E/M conservation: ResidualNucleus is added (ResNuc)
 // Class Description
 // Final state production model for theoretical models of hadron inelastic
@@ -42,5 +43,5 @@
 
 
-#include "G4HadTmpUtil.hh"		//lrint
+#include "G4HadTmpUtil.hh"  //lrint
 
 //#define debug_scatter
@@ -48,5 +49,5 @@
 G4QuasiElasticChannel::G4QuasiElasticChannel()
 {
-	theQuasiElastic=G4QuasiFreeRatios::GetPointer();
+ theQuasiElastic=G4QuasiFreeRatios::GetPointer();
 }
 
@@ -55,83 +56,129 @@
 
 G4double G4QuasiElasticChannel::GetFraction(G4Nucleus &theNucleus,
-				const G4DynamicParticle & thePrimary)
+    const G4DynamicParticle & thePrimary)
 {
-	std::pair<G4double,G4double> ratios;
-	ratios=theQuasiElastic->GetRatios(thePrimary.GetTotalMomentum(),
-			thePrimary.GetDefinition()->GetPDGEncoding(),
-			G4lrint(theNucleus.GetZ()),
-			G4lrint(theNucleus.GetN()-theNucleus.GetZ()));
-#ifdef debug_getFraction			
-	G4cout << "G4QuasiElasticChannel::ratios " << ratios.first << " x " <<ratios.second
-	       << "  = " << ratios.first*ratios.second << G4endl;
+  std::pair<G4double,G4double> ratios;
+#ifdef debug_scatter   
+  G4cout << "G4QuasiElasticChannel:: P=" << thePrimary.GetTotalMomentum()
+         << ", pPDG=" << thePrimary.GetDefinition()->GetPDGEncoding()
+         << ", Z = "  << G4lrint(theNucleus.GetZ())
+         << ", N = "  << G4lrint(theNucleus.GetN()-theNucleus.GetZ()) << G4endl;
 #endif
-	       
-	return ratios.first*ratios.second;       		
+  ratios=theQuasiElastic->GetRatios(thePrimary.GetTotalMomentum(),
+                                    thePrimary.GetDefinition()->GetPDGEncoding(),
+                                    G4lrint(theNucleus.GetZ()),
+                                    G4lrint(theNucleus.GetN()-theNucleus.GetZ()));
+#ifdef debug_scatter   
+  G4cout << "G4QuasiElasticChannel::ratios " << ratios.first << " x " <<ratios.second
+         << "  = " << ratios.first*ratios.second << G4endl;
+#endif
+        
+  return ratios.first*ratios.second;
+  //return 0.; // Switch off quasielastic (temporary) M.K.
+  //return 1.; // Only quasielastic (temporary) M.K. (DANGEROSE! Crashes at A=1!)
 }
 
 G4KineticTrackVector * G4QuasiElasticChannel::Scatter(G4Nucleus &theNucleus,
-					const G4DynamicParticle & thePrimary)
+                                                      const G4DynamicParticle & thePrimary)
 {
-	
-	
-	G4int A=G4lrint(theNucleus.GetN());
-//	G4int Z=G4lrint(theNucleus.GetZ());
-//   build Nucleus and choose random nucleon to scatter with
-	the3DNucleus.Init(theNucleus.GetN(),theNucleus.GetZ());
-	const std::vector<G4Nucleon *> nucleons=the3DNucleus.GetNucleons();
-	
-        G4int index;
-        do {
-	   index=G4lrint((A-1)*G4UniformRand());
-	   
-	} while (index < 0 || index >= static_cast<G4int>(nucleons.size()));
-	    
-//	G4double an=G4UniformRand()*A;
-	G4ParticleDefinition * pDef= nucleons[index]->GetDefinition();
-	
+  G4int A=G4lrint(theNucleus.GetN());
+  G4int Z=G4lrint(theNucleus.GetZ());                                 // M.K. ResNuc
+  //   build Nucleus and choose random nucleon to scatter with
+  the3DNucleus.Init(theNucleus.GetN(),theNucleus.GetZ());
 #ifdef debug_scatter
-	G4cout << " neutron - proton? A, Z, an, pdg" <<" "
-	       << A <<" "<<G4lrint(theNucleus.GetZ())
-	       << " "<<an <<" " << pDef->GetParticleName()<< G4endl;
+  G4cout<<"G4QElC::Scat: Before GetNucleons " << G4endl;
 #endif
-//	G4LorentzVector pNucleon(G4ThreeVector(0,0,0),pDef->GetPDGMass());
-	G4LorentzVector pNucleon=nucleons[index]->Get4Momentum();
-	
-	std::pair<G4LorentzVector,G4LorentzVector> result;
-	
-	result=theQuasiElastic->Scatter(pDef->GetPDGEncoding(),pNucleon,
-	                       thePrimary.GetDefinition()->GetPDGEncoding(),
-			        thePrimary.Get4Momentum());
-				
-	if (result.first.e() == 0.)
-	{
-	   G4cout << "Warning - G4QuasiElasticChannel::Scatter no scattering" << G4endl;
-	   return 0;       //no scatter
-	}
+  const std::vector<G4Nucleon *> nucleons=the3DNucleus.GetNucleons();
+  G4double targetNucleusMass=the3DNucleus.GetMass();                  // M.K. ResNuc
+#ifdef debug_scatter
+  G4cout<<"G4QElC::Scat: targetMass = " << targetNucleusMass << G4endl;
+#endif
+  G4LorentzVector targetNucleus4Mom(0.,0.,0.,targetNucleusMass);      // M.K. ResNuc
+  G4int index;
+  do
+  {
+    index=G4lrint((A-1)*G4UniformRand());
+  } while (index < 0 || index >= static_cast<G4int>(nucleons.size()));
+#ifdef debug_scatter
+  G4cout<<"G4QElC::Scat: index = " << index << G4endl;
+#endif
+  G4ParticleDefinition * pDef= nucleons[index]->GetDefinition();
+  G4int resA=A-1;                                                     // M.K. ResNuc
+  G4int resZ=Z-static_cast<int>(pDef->GetPDGCharge());                // M.K. ResNuc
+  G4ParticleDefinition* resDef=G4Neutron::Neutron(); // Resolve t-p=nn problem M.K. ResNuc
+  G4double residualNucleusMass=resDef->GetPDGMass();                  // M.K. ResNuc
+  if(resZ)
+  {
+    resDef=G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ);// M.K. ResNuc
+    residualNucleusMass=resDef->GetPDGMass();                         // M.K. ResNuc
+  }
+  else residualNucleusMass*=resA;                           // resA=resN M.K. ResNuc 
+#ifdef debug_scatter
+  G4cout<<"G4QElChan::Scatter: neutron - proton? A ="<<A<<", Z="<<Z<<", projName="
+        <<pDef->GetParticleName()<<G4endl;
+#endif
+  // G4LorentzVector pNucleon(G4ThreeVector(0,0,0),pDef->GetPDGMass());
+  G4LorentzVector pNucleon=nucleons[index]->Get4Momentum();
+  G4double residualNucleusEnergy=std::sqrt(residualNucleusMass*residualNucleusMass+
+                                           pNucleon.vect().mag2());   // M.K. ResNuc
+  pNucleon.setE(targetNucleusMass-residualNucleusEnergy);             // M.K. ResNuc
+  G4LorentzVector residualNucleus4Mom=targetNucleus4Mom-pNucleon;     // M.K. ResNuc
+ 
+  std::pair<G4LorentzVector,G4LorentzVector> result;
+ 
+  result=theQuasiElastic->Scatter(pDef->GetPDGEncoding(),pNucleon,
+                                  thePrimary.GetDefinition()->GetPDGEncoding(),
+  thePrimary.Get4Momentum());
+  G4LorentzVector scatteredHadron4Mom=result.second;                  // M.K. ResNuc
+  if (result.first.e() <= 0.)
+  {
+    //G4cout << "Warning - G4QuasiElasticChannel::Scatter no scattering" << G4endl;
+    //return 0;       //no scatter
+    G4LorentzVector scatteredHadron4Mom=thePrimary.Get4Momentum();   // M.K. ResNuc
+    residualNucleus4Mom=G4LorentzVector(0.,0.,0.,targetNucleusMass); // M.K. ResNuc
+    resDef=G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);    // M.K. ResNuc
+  }
 
 #ifdef debug_scatter
-        G4LorentzVector EpConservation=pNucleon+thePrimary.Get4Momentum() 
-	                                - result.first - result.second;
-	if (   (EpConservation.vect().mag2() > 0.01*MeV*MeV )
-	    || 	(std::abs(EpConservation.e()) > 0.1 * MeV ) ) 
-	{
-	    G4cout << "Warning - G4QuasiElasticChannel::Scatter E-p non conservation : "
-	            << EpConservation << G4endl;
-	}    
+  G4LorentzVector EpConservation=pNucleon+thePrimary.Get4Momentum() 
+                                 - result.first - result.second;
+  if (   (EpConservation.vect().mag2() > .01*MeV*MeV )
+      || (std::abs(EpConservation.e()) > 0.1 * MeV ) ) 
+  {
+    G4cout << "Warning - G4QuasiElasticChannel::Scatter E-p non conservation : "
+           << EpConservation << G4endl;
+  }    
 #endif
 
-	G4KineticTrackVector * ktv;
-	ktv=new G4KineticTrackVector();
-	G4KineticTrack * sPrim=new G4KineticTrack(thePrimary.GetDefinition(),
-			0.,G4ThreeVector(0), result.second);
-	ktv->push_back(sPrim);
-	G4KineticTrack * sNuc=new G4KineticTrack(pDef,
-			0.,G4ThreeVector(0), result.first);
-	ktv->push_back(sNuc);
-	
+  G4KineticTrackVector * ktv;
+  ktv=new G4KineticTrackVector();
+  G4KineticTrack * sPrim=new G4KineticTrack(thePrimary.GetDefinition(),
+                                            0.,G4ThreeVector(0), scatteredHadron4Mom);
+  ktv->push_back(sPrim);
+  if (result.first.e() > 0.)
+  {
+    G4KineticTrack * sNuc=new G4KineticTrack(pDef, 0.,G4ThreeVector(0), result.first);
+    ktv->push_back(sNuc);
+  }
+  if(resZ || resA==1) // For the only neutron or for tnuclei with Z>0    M.K. ResNuc
+  {
+    G4KineticTrack * rNuc=new G4KineticTrack(resDef,
+                           0.,G4ThreeVector(0), residualNucleus4Mom); // M.K. ResNuc
+    ktv->push_back(rNuc);                                             // M.K. ResNuc
+  }
+  else // The residual nucleus consists of only neutrons                 M.K. ResNuc
+  {
+    residualNucleus4Mom/=resA;     // Split 4-mom of A*n system equally  M.K. ResNuc
+    for(G4int in=0; in<resA; in++) // Loop over neutrons in A*n system.  M.K. ResNuc
+    {
+      G4KineticTrack* rNuc=new G4KineticTrack(resDef,
+                           0.,G4ThreeVector(0), residualNucleus4Mom); // M.K. ResNuc
+      ktv->push_back(rNuc);                                           // M.K. ResNuc
+    }
+  }
 #ifdef debug_scatter
-	G4cout << " scattered Nucleon : " << result.first << " mass " <<result.first.mag() << G4endl;
-	G4cout << " scattered Project : " << result.second << " mass " <<result.second.mag() << G4endl;
+  G4cout<<"G4QElC::Scat: Nucleon: "<<result.first <<" mass "<<result.first.mag() << G4endl;
+  G4cout<<"G4QElC::Scat: Project: "<<result.second<<" mass "<<result.second.mag()<< G4endl;
 #endif
-	return ktv;			       
+  return ktv;
 }
Index: trunk/source/processes/hadronic/models/theo_high_energy/src/G4TheoFSGenerator.cc
===================================================================
--- trunk/source/processes/hadronic/models/theo_high_energy/src/G4TheoFSGenerator.cc	(revision 1007)
+++ trunk/source/processes/hadronic/models/theo_high_energy/src/G4TheoFSGenerator.cc	(revision 1055)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4TheoFSGenerator.cc,v 1.9 2007/11/13 16:01:36 gunter Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4TheoFSGenerator.cc,v 1.11 2009/04/09 08:28:42 mkossov Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // G4TheoFSGenerator
@@ -91,5 +91,40 @@
      if ( theQuasielastic->GetFraction(theNucleus, *aPart) > G4UniformRand() )
      {
-        G4KineticTrackVector *result= theQuasielastic->Scatter(theNucleus, *aPart);
+       //G4cout<<"___G4TheoFSGenerator: before Scatter (1) QE=" << theQuasielastic<<G4endl;
+       G4KineticTrackVector *result= theQuasielastic->Scatter(theNucleus, *aPart);
+       //G4cout << "^^G4TheoFSGenerator: after Scatter (1) " << G4endl;
+       if (result)
+       {
+	    for(unsigned int  i=0; i<result->size(); i++)
+	    {
+	      G4DynamicParticle * aNew = 
+		 new G4DynamicParticle(result->operator[](i)->GetDefinition(),
+                        	       result->operator[](i)->Get4Momentum().e(),
+                        	       result->operator[](i)->Get4Momentum().vect());
+	      theParticleChange->AddSecondary(aNew);
+	      delete result->operator[](i);
+	    }
+	    delete result;
+	   
+       } else 
+       {
+	    theParticleChange->SetStatusChange(isAlive);
+	    theParticleChange->SetEnergyChange(thePrimary.GetKineticEnergy());
+	    theParticleChange->SetMomentumChange(thePrimary.Get4Momentum().vect().unit());
+ 
+       }
+	return theParticleChange;
+     } 
+  }
+  if ( theProjectileDiffraction) {
+  
+     if ( theProjectileDiffraction->GetFraction(theNucleus, *aPart) > G4UniformRand() )
+      // strictly, this returns fraction on inelastic, so quasielastic should 
+	//    already be substracted, ie. quasielastic should always be used
+	//     before diffractive
+     {
+       //G4cout << "____G4TheoFSGenerator: before Scatter (2) " << G4endl;
+       G4KineticTrackVector *result= theProjectileDiffraction->Scatter(theNucleus, *aPart);
+       //G4cout << "^^^^G4TheoFSGenerator: after Scatter (2) " << G4endl;
 	if (result)
 	{
@@ -115,39 +150,8 @@
      } 
   }
-  if ( theProjectileDiffraction) {
-  
-     if ( theProjectileDiffraction->GetFraction(theNucleus, *aPart) > G4UniformRand() )
-        // strictly, this returns fraction on inelastic, so quasielastic should 
-	//    already be substracted, ie. quasielastic should always be used
-	//     before diffractive
-     {
-        G4KineticTrackVector *result= theProjectileDiffraction->Scatter(theNucleus, *aPart);
-	if (result)
-	{
-	    for(unsigned int  i=0; i<result->size(); i++)
-	    {
-	      G4DynamicParticle * aNew = 
-		 new G4DynamicParticle(result->operator[](i)->GetDefinition(),
-                        	       result->operator[](i)->Get4Momentum().e(),
-                        	       result->operator[](i)->Get4Momentum().vect());
-	      theParticleChange->AddSecondary(aNew);
-	      delete result->operator[](i);
-	    }
-	    delete result;
-	   
-	} else 
-	{
-	    theParticleChange->SetStatusChange(isAlive);
-	    theParticleChange->SetEnergyChange(thePrimary.GetKineticEnergy());
-	    theParticleChange->SetMomentumChange(thePrimary.Get4Momentum().vect().unit());
- 
-	}
-	return theParticleChange;
-     } 
-  }
-
+  //G4cout << "____G4TheoFSGenerator: before Scatter (3) " << G4endl;
   G4KineticTrackVector * theInitialResult =
                theHighEnergyGenerator->Scatter(theNucleus, *aPart);
-  
+  //G4cout << "^^^^G4TheoFSGenerator: after Scatter (3) " << G4endl;
   
   G4ReactionProductVector * theTransportResult = NULL;
