lars/geant4
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Import Geant4 9.2.0 source tree

Browse Source
This commit is contained in:
Gabriele Cosmo
2016-06-09 15:58:43 +02:00
parent 96c8bcd0af
commit b79225fb37
7544 changed files with 245407 additions and 91099 deletions
Download Patch File Download Diff File Expand all files Collapse all files
source/processes/hadronic/models/coherent_elastic/src/G4ChargeExchange.cc
+56 -92
View File
@@ -24,8 +24,8 @@
// ********************************************************************
//
//
// $Id: G4ChargeExchange.cc,v 1.11 2007/05/25 17:46:52 dennis Exp $
// GEANT4 tag $Name: geant4-09-01 $
// $Id: G4ChargeExchange.cc,v 1.14 2008/12/18 13:01:48 gunter Exp $
// GEANT4 tag $Name: geant4-09-02 $
//
//
// G4 Model: Charge and strangness exchange based on G4LightMedia model
@@ -40,32 +40,16 @@
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4IonTable.hh"
#include "G4QElasticCrossSection.hh"
#include "G4VQCrossSection.hh"
#include "G4ElasticHadrNucleusHE.hh"
#include "Randomize.hh"
#include "G4HadronElastic.hh"
#include "G4NucleiProperties.hh"
G4ChargeExchange::G4ChargeExchange(G4HadronElastic* hel, G4double elim,
G4double ehigh)
: G4HadronicInteraction("G4ChargeExchange"),
fElastic(hel),
native(false),
ekinlim(elim),
ekinhigh(ehigh)
G4ChargeExchange::G4ChargeExchange() : G4HadronicInteraction("Charge Exchange")
{
SetMinEnergy( 0.0*GeV );
SetMaxEnergy( 100.*TeV );
ekinlow = 19.0*MeV;
verboseLevel= 0;
if(!fElastic) {
native = true;
fElastic = new G4HadronElastic();
}
qCManager = fElastic->GetCS();
hElastic = fElastic->GetHElastic();
lowEnergyRecoilLimit = 100.*keV;
lowestEnergyLimit = 1.*MeV;
theProton = G4Proton::Proton();
theNeutron = G4Neutron::Neutron();
@@ -99,9 +83,7 @@ G4ChargeExchange::G4ChargeExchange(G4HadronElastic* hel, G4double elim,
}
G4ChargeExchange::~G4ChargeExchange()
{
if(native) delete fElastic;
}
{}
G4HadFinalState* G4ChargeExchange::ApplyYourself(
const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
@@ -116,7 +98,7 @@ G4HadFinalState* G4ChargeExchange::ApplyYourself(
G4int Z = static_cast<G4int>(zTarget+0.5);
G4int A = static_cast<G4int>(aTarget+0.5);
if(ekin == 0.0 || A < 3) {
if(ekin <= lowestEnergyLimit || A < 3) {
theParticleChange.SetEnergyChange(ekin);
theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
return &theParticleChange;
@@ -132,7 +114,6 @@ G4HadFinalState* G4ChargeExchange::ApplyYourself(
// Scattered particle referred to axis of incident particle
const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
G4double m1 = theParticle->GetPDGMass();
G4int N = A - Z;
G4int projPDG = theParticle->GetPDGEncoding();
@@ -144,12 +125,7 @@ G4HadFinalState* G4ChargeExchange::ApplyYourself(
G4ParticleDefinition * theDef = 0;
if (Z == 1 && A == 3) theDef = theT;
else if (Z == 2 && A == 3) theDef = theHe3;
else if (Z == 2 && A == 4) theDef = theA;
else theDef = G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);
G4double m2 = theDef->GetPDGMass();
G4double m2 = G4NucleiProperties::GetNuclearMass((G4double)A, (G4double)Z);
G4LorentzVector lv1 = aParticle->Get4Momentum();
G4LorentzVector lv0(0.0,0.0,0.0,m2);
@@ -258,80 +234,36 @@ G4HadFinalState* G4ChargeExchange::ApplyYourself(
else theRecoil = theDef;
G4double etot = lv0.e() + lv1.e();
if(etot < m11 + m21) return &theParticleChange;
// kinematiacally impossible
if(etot < m11 + m21) {
theParticleChange.SetEnergyChange(ekin);
theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
return &theParticleChange;
}
G4ThreeVector p1 = lv1.vect();
G4double e1 = 0.5*etot*(1.0 + (m21*m21 - m11*m11)/(etot*etot));
G4double e1 = 0.5*etot*(1.0 - (m21*m21 - m11*m11)/(etot*etot));
// G4double e2 = etot - e1;
G4double ptot = std::sqrt(e1*e1 - m11*m11);
G4double tmax = 4.0*ptot*ptot;
G4double t = 0.0;
G4double g2 = GeV*GeV;
// Choose generator
G4ElasticGenerator gtype = fLElastic;
if ((theParticle == theProton || theParticle == theNeutron) &&
Z <= 2 && ekin >= ekinlow) {
gtype = fQElastic;
} else {
if(ekin >= ekinlow) gtype = fSWave;
else if(ekin >= ekinhigh) gtype = fHElastic;
}
// Sample t
if(gtype == fQElastic) {
if (verboseLevel > 1)
G4cout << "G4ChargeExchange: Z= " << Z << " N= "
<< N << " pdg= " << projPDG
<< " mom(GeV)= " << plab/GeV << " " << qCManager << G4endl;
if(Z == 1 && N == 2) N = 1;
else if (Z == 2 && N == 1) N = 2;
G4double cs = qCManager->GetCrossSection(false,plab,Z,N,projPDG);
if(cs > 0.0) t = qCManager->GetExchangeT(Z,N,projPDG);
else gtype = fSWave;
}
if(gtype == fHElastic) {
t = hElastic->SampleT(theParticle,plab,Z,A);
if(t > tmax) gtype = fSWave;
}
if(gtype == fLElastic) {
t = GeV*GeV*fElastic->SampleT(ptot,m1,m2,aTarget);
if(t > tmax) gtype = fSWave;
}
// NaN finder
if(!(t < 0.0 || t >= 0.0)) {
if (verboseLevel > -1) {
G4cout << "G4ChargeExchange:WARNING: Z= " << Z << " N= "
<< N << " pdg= " << projPDG
<< " mom(GeV)= " << plab/GeV
<< " the model type " << gtype;
if(gtype == fQElastic) G4cout << " CHIPS ";
else if(gtype == fLElastic) G4cout << " LElastic ";
else if(gtype == fHElastic) G4cout << " HElastic ";
G4cout << " t= " << t
<< " S-wave will be sampled"
<< G4endl;
}
gtype = fSWave;
}
if(gtype == fSWave) t = G4UniformRand()*tmax;
G4double t = g2*SampleT(tmax/g2,aTarget);
if(verboseLevel>1)
G4cout <<"type= " << gtype <<" t= " << t << " tmax= " << tmax
G4cout <<"## G4ChargeExchange t= " << t << " tmax= " << tmax
<< " ptot= " << ptot << G4endl;
// Sampling in CM system
G4double phi = G4UniformRand()*twopi;
G4double cost = 1. - 2.0*t/tmax;
if(std::abs(cost) > 1.0) cost = -1.0 + 2.0*G4UniformRand();
if(std::abs(cost) > 1.0) cost = 1.0;
G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
if (verboseLevel > 1)
G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
//if (verboseLevel > 1)
// G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
v1 *= ptot;
@@ -342,20 +274,52 @@ G4HadFinalState* G4ChargeExchange::ApplyYourself(
nlv1.boost(bst);
theParticleChange.SetStatusChange(stopAndKill);
theParticleChange.SetEnergyChange(0.0);
G4DynamicParticle * aSec = new G4DynamicParticle(theSecondary, nlv1);
theParticleChange.AddSecondary(aSec);
G4double erec = nlv0.e() - m21;
//G4cout << "erec= " <<erec << " Esec= " << aSec->GetKineticEnergy() << G4endl;
if(theHyperon) {
theParticleChange.SetLocalEnergyDeposit(erec);
aSec = new G4DynamicParticle();
aSec->SetDefinition(theRecoil);
aSec->SetKineticEnergy(0.0);
} else if(erec > ekinlim) {
} else if(erec > lowEnergyRecoilLimit) {
aSec = new G4DynamicParticle(theRecoil, nlv0);
theParticleChange.AddSecondary(aSec);
} else {
if(erec < 0.0) erec = 0.0;
theParticleChange.SetLocalEnergyDeposit(erec);
}
return &theParticleChange;
}
G4double G4ChargeExchange::SampleT(G4double tmax, G4double A)
{
G4double aa, bb, cc, dd;
if (A <= 62.) {
aa = std::pow(A, 1.63);
bb = 14.5*std::pow(A, 0.66);
cc = 1.4*std::pow(A, 0.33);
dd = 10.;
} else {
aa = std::pow(A, 1.33);
bb = 60.*std::pow(A, 0.33);
cc = 0.4*std::pow(A, 0.40);
dd = 10.;
}
G4double x1 = (1.0 - std::exp(-tmax*bb))*aa/bb;
G4double x2 = (1.0 - std::exp(-tmax*dd))*cc/dd;
G4double t;
G4double y = bb;
if(G4UniformRand()*(x1 + x2) < x2) y = dd;
do {t = -std::log(G4UniformRand())/y;} while (t > tmax);
return t;
}
source/processes/hadronic/models/coherent_elastic/src/G4ChargeExchangeProcess.cc
+23 -163
View File
@@ -24,11 +24,11 @@
// ********************************************************************
//
//
// $Id: G4ChargeExchangeProcess.cc,v 1.9 2007/01/30 10:23:26 vnivanch Exp $
// GEANT4 tag $Name: geant4-09-01 $
// $Id: G4ChargeExchangeProcess.cc,v 1.15 2008/11/27 16:43:00 vnivanch Exp $
// GEANT4 tag $Name: geant4-09-02 $
//
//
// Geant4 Hadron Elastic Scattering Process -- header file
// Geant4 Hadron Charge Exchange Process -- source file
//
// Created 21 April 2006 V.Ivanchenko
//
@@ -44,23 +44,19 @@
#include "globals.hh"
#include "G4CrossSectionDataStore.hh"
#include "G4HadronElasticDataSet.hh"
#include "G4VQCrossSection.hh"
#include "G4QElasticCrossSection.hh"
#include "G4QCHIPSWorld.hh"
#include "G4Element.hh"
#include "G4ElementVector.hh"
#include "G4IsotopeVector.hh"
#include "G4Neutron.hh"
#include "G4Proton.hh"
#include "G4HadronElastic.hh"
#include "G4PhysicsLinearVector.hh"
G4ChargeExchangeProcess::G4ChargeExchangeProcess(const G4String& procName)
: G4HadronicProcess(procName), first(true)
{
thEnergy = 19.*MeV;
SetProcessSubType(fChargeExchange);
thEnergy = 20.*MeV;
verboseLevel= 1;
qCManager = 0;
AddDataSet(new G4HadronElasticDataSet);
theProton = G4Proton::Proton();
theNeutron = G4Neutron::Neutron();
@@ -98,11 +94,6 @@ G4ChargeExchangeProcess::~G4ChargeExchangeProcess()
delete factors;
}
void G4ChargeExchangeProcess::SetQElasticCrossSection(G4VQCrossSection* p)
{
qCManager = p;
}
void G4ChargeExchangeProcess::
BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
{
@@ -119,63 +110,23 @@ BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
theParticle == theK0S || theParticle == theK0L) {
G4double F[n] = {0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.1,0.09,0.07};
factors = new G4PhysicsLinearVector(0.0,1.8*GeV,n);
factors = new G4PhysicsLinearVector(0.0,2.0*GeV,n);
for(size_t i=0; i<n; i++) {factors->PutValue(i,F[i]);}
} else {
G4double F[n] = {0.50,0.45,0.40,0.35,0.30,0.25,0.06,0.04,0.005,0.0};
factors = new G4PhysicsLinearVector(0.0,3.6*GeV,n);
factors = new G4PhysicsLinearVector(0.0,4.0*GeV,n);
for(size_t i=0; i<n; i++) {factors->PutValue(i,F[i]);}
}
//factors->SetSpline(true);
if(verboseLevel>1)
G4cout << "G4ChargeExchangeProcess for "
<< theParticle->GetParticleName()
<< G4endl;
}
store->BuildPhysicsTable(aParticleType);
}
G4double G4ChargeExchangeProcess::GetMeanFreePath(const G4Track& track,
G4double,
G4ForceCondition* cond)
{
*cond = NotForced;
const G4DynamicParticle* dp = track.GetDynamicParticle();
const G4Material* material = track.GetMaterial();
cross = 0.0;
G4double x = DBL_MAX;
// The process is effective only above the threshold
if(dp->GetKineticEnergy() < thEnergy) return x;
// Compute cross sesctions
const G4ElementVector* theElementVector = material->GetElementVector();
const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
G4double temp = material->GetTemperature();
G4int nelm = material->GetNumberOfElements();
if(verboseLevel>1)
G4cout << "G4ChargeExchangeProcess get mfp for "
<< theParticle->GetParticleName()
<< " p(GeV)= " << dp->GetTotalMomentum()/GeV
<< " in " << material->GetName()
<< G4endl;
for (G4int i=0; i<nelm; i++) {
const G4Element* elm = (*theElementVector)[i];
G4double x = GetMicroscopicCrossSection(dp, elm, temp);
cross += theAtomNumDensityVector[i]*x;
xsec[i] = cross;
}
if(verboseLevel>1)
G4cout << "G4ChargeExchangeProcess cross(1/mm)= " << cross
<< " E(MeV)= " << dp->GetKineticEnergy()
<< " " << theParticle->GetParticleName()
<< " in " << material->GetName()
<< G4endl;
if(cross > DBL_MIN) x = 1./cross;
return x;
G4HadronicProcess::BuildPhysicsTable(aParticleType);
}
G4double G4ChargeExchangeProcess::GetMicroscopicCrossSection(
@@ -187,7 +138,9 @@ G4double G4ChargeExchangeProcess::GetMicroscopicCrossSection(
G4double Z = elm->GetZ();
G4int iz = G4int(Z);
G4double x = 0.0;
if(iz == 1) return x;
// The process is effective only above the threshold
if(iz == 1 || dp->GetKineticEnergy() < thEnergy) return x;
if(verboseLevel>1)
G4cout << "G4ChargeExchangeProcess compute GHAD CS for element "
@@ -195,19 +148,6 @@ G4double G4ChargeExchangeProcess::GetMicroscopicCrossSection(
<< G4endl;
x = store->GetCrossSection(dp, elm, temp);
// NaN finder
if(!(x < 0.0 || x >= 0.0)) {
if (verboseLevel > -1) {
G4cout << "G4ChargeExchangeProcess WARNING: Z= " << iz
<< " pdg= " << pPDG
<< " mom(GeV)= " << dp->GetTotalMomentum()/GeV
<< " cross= " << x
<< " set to zero"
<< G4endl;
}
x = 0.0;
}
if(verboseLevel>1)
G4cout << "G4ChargeExchangeProcess cross(mb)= " << x/millibarn
<< " E(MeV)= " << dp->GetKineticEnergy()
@@ -216,106 +156,26 @@ G4double G4ChargeExchangeProcess::GetMicroscopicCrossSection(
<< G4endl;
G4bool b;
G4double A = elm->GetN();
x *= factors->GetValue(dp->GetTotalMomentum(), b)/std::pow(A, 0.42);
G4double ptot = dp->GetTotalMomentum();
x *= factors->GetValue(ptot, b)/std::pow(A, 0.42);
if(theParticle == thePiPlus || theParticle == theProton ||
theParticle == theKPlus || theParticle == theANeutron)
x *= (1.0 - Z/A);
{ x *= (1.0 - Z/A); }
else if(theParticle == thePiMinus || theParticle == theNeutron ||
theParticle == theKMinus || theParticle == theAProton)
x *= Z/A;
{ x *= Z/A; }
if(theParticle->GetPDGMass() < GeV) {
if(ptot > 2.*GeV) x *= 4.0*GeV*GeV/(ptot*ptot);
}
if(verboseLevel>1)
G4cout << "Corrected cross(mb)= " << x/millibarn << G4endl;
return x;
}
G4VParticleChange* G4ChargeExchangeProcess::PostStepDoIt(
const G4Track& track,
const G4Step& step)
{
G4ForceCondition* cn = 0;
aParticleChange.Initialize(track);
G4double mfp = GetMeanFreePath(track, 0.0, cn);
if(mfp == DBL_MAX) return G4VDiscreteProcess::PostStepDoIt(track,step);
G4double kineticEnergy = track.GetKineticEnergy();
G4Material* material = track.GetMaterial();
// Select element
const G4ElementVector* theElementVector = material->GetElementVector();
G4Element* elm = (*theElementVector)[0];
G4int nelm = material->GetNumberOfElements() - 1;
if (nelm > 0) {
G4double x = G4UniformRand()*cross;
G4int i = -1;
do {i++;} while (x > xsec[i] && i < nelm);
elm = (*theElementVector)[i];
}
G4double Z = elm->GetZ();
G4double A = G4double(G4int(elm->GetN()+0.5));
// Select isotope
G4IsotopeVector* isv = elm->GetIsotopeVector();
G4int ni = 0;
if(isv) ni = isv->size();
if(ni == 1) {
A = G4double((*isv)[0]->GetN());
} else if(ni > 1) {
G4double* ab = elm->GetRelativeAbundanceVector();
G4int j = -1;
ni--;
G4double y = G4UniformRand();
do {
j++;
y -= ab[j];
} while (y > 0.0 && j < ni);
A = G4double((*isv)[j]->GetN());
}
G4HadronicInteraction* hadi =
ChooseHadronicInteraction( kineticEnergy, material, elm);
// Initialize the hadronic projectile from the track
G4HadProjectile thePro(track);
if(verboseLevel>1)
G4cout << "G4ChargeExchangeProcess::PostStepDoIt for "
<< theParticle->GetParticleName()
<< " Target Z= " << Z
<< " A= " << A << G4endl;
targetNucleus.SetParameters(A, Z);
aParticleChange.Initialize(track);
G4HadFinalState* result = hadi->ApplyYourself(thePro, targetNucleus);
G4ThreeVector indir = track.GetMomentumDirection();
G4int nsec = result->GetNumberOfSecondaries();
if(verboseLevel>1)
G4cout << "Efin= " << result->GetEnergyChange()
<< " de= " << result->GetLocalEnergyDeposit()
<< " nsec= " << nsec
<< G4endl;
if(nsec > 0) {
aParticleChange.ProposeEnergy(0.0);
aParticleChange.ProposeTrackStatus(fStopAndKill);
aParticleChange.ProposeLocalEnergyDeposit(result->GetLocalEnergyDeposit());
aParticleChange.SetNumberOfSecondaries(nsec);
for(G4int j=0; j<nsec; j++) {
G4DynamicParticle* p = result->GetSecondary(j)->GetParticle();
G4ThreeVector pdir = p->GetMomentumDirection();
// G4cout << "recoil " << pdir << G4endl;
pdir = pdir.rotateUz(indir);
// G4cout << "recoil rotated " << pdir << G4endl;
p->SetMomentumDirection(pdir);
aParticleChange.AddSecondary(p);
}
}
result->Clear();
return G4VDiscreteProcess::PostStepDoIt(track,step);
}
G4bool G4ChargeExchangeProcess::
IsApplicable(const G4ParticleDefinition& aParticleType)
{
source/processes/hadronic/models/coherent_elastic/src/G4DiffuseElastic.cc
+2 -5
View File
@@ -23,8 +23,8 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4DiffuseElastic.cc,v 1.18 2007/11/14 18:42:29 vnivanch Exp $
// GEANT4 tag $Name: geant4-09-01 $
// $Id: G4DiffuseElastic.cc,v 1.20 2008/01/14 10:39:13 vnivanch Exp $
// GEANT4 tag $Name: geant4-09-02 $
//
//
// Physics model class G4DiffuseElastic
@@ -39,9 +39,6 @@
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4IonTable.hh"
#include "G4QElasticCrossSection.hh"
#include "G4VQCrossSection.hh"
#include "G4ElasticHadrNucleusHE.hh"
#include "Randomize.hh"
#include "G4Integrator.hh"
source/processes/hadronic/models/coherent_elastic/src/G4ElasticHadrNucleusHE.cc
+201 -358
View File
@@ -24,8 +24,8 @@
// ********************************************************************
//
//
// $Id: G4ElasticHadrNucleusHE.cc,v 1.75 2007/11/13 17:22:51 vnivanch Exp $
// GEANT4 tag $Name: geant4-09-01 $
// $Id: G4ElasticHadrNucleusHE.cc,v 1.79 2008/01/14 10:39:13 vnivanch Exp $
// GEANT4 tag $Name: geant4-09-02 $
//
//
// The generator of high energy hadron-nucleus elastic scattering
@@ -400,9 +400,8 @@ G4HadFinalState * G4ElasticHadrNucleusHE::ApplyYourself(
sint = std::sqrt((1.0-cost)*(1.0+cost));
}
if (verboseLevel>1)
{
G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
}
G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
v1 *= ptot;
G4LorentzVector nlv1( v1.x(), v1.y(), v1.z(), std::sqrt(ptot*ptot + m1*m1));
@@ -468,7 +467,7 @@ G4double G4ElasticHadrNucleusHE::
{
G4double plab = inLabMom/GeV; // (GeV/c)
G4double Q2 = 0;
// G4int ii, jj;
iHadrCode = p->GetPDGEncoding();
NumbN = N;
@@ -502,54 +501,46 @@ G4double G4ElasticHadrNucleusHE::
G4double T = sqrt(plab*plab+hMass2)-hMass;
if(T < 0.4) Q2 = 0;
else Q2 = HadronProtonQ2(p, plab);
if(T > 0.4) Q2 = HadronProtonQ2(p, plab);
if (verboseLevel>1)
G4cout<<" Proton : Q2 "<<Q2<<G4endl;
return Q2*GeV2;
if (verboseLevel>1)
G4cout<<" Proton : Q2 "<<Q2<<G4endl;
}
else
{
G4ElasticData* ElD1 = SetOfElasticData[idx][Z];
// Construct elastic data
if(!ElD1)
{
G4double AWeight = nistManager->GetAtomicMassAmu(Z);
ElD1 = new G4ElasticData(p, Z, AWeight, Energy);
SetOfElasticData[idx][Z] = ElD1;
if(verboseLevel > 1)
{
G4cout<< " G4ElasticHadrNucleusHE::SampleT: new record " << idx
<< " for " << p->GetParticleName() << " Z= " << Z
<< G4endl;
G4ElasticData* ElD1 = SetOfElasticData[idx][Z];
// Construct elastic data
if(!ElD1)
{
G4double AWeight = nistManager->GetAtomicMassAmu(Z);
ElD1 = new G4ElasticData(p, Z, AWeight, Energy);
SetOfElasticData[idx][Z] = ElD1;
if(verboseLevel > 1)
{
G4cout<< " G4ElasticHadrNucleusHE::SampleT: new record " << idx
<< " for " << p->GetParticleName() << " Z= " << Z
<< G4endl;
}
}
hMass = ElD1->massGeV;
hMass2 = ElD1->mass2GeV2;
G4double M = ElD1->massA;
G4double M2 = ElD1->massA2;
G4double plab2 = plab*plab;
G4double Q2max = 4.*plab2*M2/
(hMass2 + M2 + 2.*M*std::sqrt(plab2 + hMass2));
// sample scattering
G4double T = sqrt(plab2+hMass2)-hMass;
if(T > 0.4) Q2 = HadronNucleusQ2_2(ElD1, Z, plab, Q2max);
if(verboseLevel > 1)
G4cout<<" SampleT: Q2(GeV^2)= "<<Q2<< " t/tmax= " << Q2/Q2max <<G4endl;
}
}
hMass = ElD1->massGeV;
hMass2 = ElD1->mass2GeV2;
G4double M = ElD1->massA;
G4double M2 = ElD1->massA2;
G4double plab2 = plab*plab;
G4double Q2max = 4.*plab2*M2/(hMass2 + M2 + 2.*M*std::sqrt(plab2 + hMass2));
// sample scattering
G4double T = sqrt(plab2+hMass2)-hMass;
if(T < 0.4) Q2 = 0;
else Q2 = HadronNucleusQ2_2(ElD1, Z, plab, Q2max);
if(verboseLevel > 1)
G4cout<<" SampleT: Q2(GeV^2)= "<<Q2<< " t/tmax= " << Q2/Q2max <<G4endl;
return Q2*GeV2;
}
}
//////////////////////////////////////////////////////////////////////////
@@ -560,7 +551,6 @@ G4double G4ElasticHadrNucleusHE::
HadronNucleusQ2_2(G4ElasticData* pElD, G4int Z,
G4double plab, G4double tmax)
{
G4double LineFq2[ONQ2];
G4double Rand = G4UniformRand();
@@ -572,11 +562,10 @@ G4double G4ElasticHadrNucleusHE::
G4double ekin = std::sqrt(hMass2 + ptot2) - hMass;
if(verboseLevel > 1)
G4cout<<"Q2_2: ekin plab "<<ekin<<" "<<plab<<" tmax "<<tmax<<G4endl;
G4cout<<"Q2_2: ekin plab "<<ekin<<" "<<plab<<" tmax "<<tmax<<G4endl;
// Find closest energy bin
G4int NumbOnE;
for( NumbOnE = 0; NumbOnE < NENERGY-1; NumbOnE++ )
{
if( ekin <= LowEdgeEnergy[NumbOnE+1] ) break;
@@ -585,156 +574,65 @@ G4double G4ElasticHadrNucleusHE::
G4int index = NumbOnE*ONQ2;
G4double Weight= 1.0;
G4double rmax = 1.0;
// Select kinematics for node energy
G4double T = Energy[NumbOnE];
hLabMomentum2 = T*(T + 2.*hMass);
G4double Q2max = pElD->maxQ2[NumbOnE];
G4int length = pElD->dnkE[NumbOnE];
G4bool isIni = false;
// Build first part of the vector
// Build vector
if(length == 0)
{
isIni = true;
R1 = pElD->R1;
R2 = pElD->R2;
Aeff = pElD->Aeff;
Pnucl = pElD->Pnucl;
hLabMomentum = std::sqrt(hLabMomentum2);
DefineHadronValues(Z);
if(verboseLevel >0)
{
G4cout<<"1 plab T "<<plab<<" "<<T<<" sigTot B ReIm "
<<HadrTot<<" "<<HadrSlope<<" "<<HadrReIm<<G4endl;
G4cout<<" R1 R2 Aeff p "<<R1<<" "<<R2<<" "<<Aeff<<" "
<<Pnucl<<G4endl;
}
// G4int AWeight = pElD->AtomicWeight;
/// Weight = GetLightFq2(Z, AWeight, Q2max);
pElD->CrossSecMaxQ2[NumbOnE] = Weight;
if(verboseLevel > 1)
G4cout<<" HadrNucleusQ2_2: NumbOnE= " << NumbOnE
<< " length= " << length
<< " Weight "<<Weight
<< " Q2max= " << Q2max
<< " ekin= " << ekin <<G4endl;
pElD->TableCrossSec[index] = 0;
dQ2 = pElD->TableQ2[1]-pElD->TableQ2[0];
GetHeavyFq2(NumbN, LineFq2); // %%%%%%%%%%%%%%%%%%%%%%%%%
for(G4int ii=0; ii<ONQ2; ii++)
{
if(verboseLevel > 2)
G4cout<<" ii LineFq2 "<<ii<<" "<<LineFq2[ii]/LineFq2[ONQ2-1]
<<" dF(q2) "<<LineFq2[ii]-LineFq2[ii-1]<<G4endl;
pElD->TableCrossSec[index+ii] = LineFq2[ii]/LineFq2[ONQ2-1];
}
pElD->dnkE[NumbOnE] = ONQ2;
rmax = 1;
/*
for(G4int ii=1; ii<ONQ0; ii++)
{
Q2 = pElD->TableQ2[ii];
if(Q2 < Q2max) Buf = GetLightFq2(Z, AWeight, Q2)/Weight;
else Buf = 1.0;
pElD->TableCrossSec[index+ii] = Buf;
if(verboseLevel > 1)
G4cout<<" HadrNucleusQ2_2: ii= " << ii << " Q2= "
<<Q2 <<" p= " <<Buf<<" B*W "<<Buf*Weight<<G4endl;
CrSec = HadrNucDifferCrSec(NumbN, Q2); //$$$$$$$$$$$$$$$$$$$
if(verboseLevel > 1)
G4cout<<" ii Q2 DiffCrSec "<<ii<<" "<<Q2<<" "<<CrSec<<G4endl;
} // for ii
rmax = Buf;
length = ONQ0;
pElD->dnkE[NumbOnE] = ONQ0;
*/
}
else
{
rmax = pElD->TableCrossSec[index+length-1];
}
G4double* dNumbFQ2 = &(pElD->TableCrossSec[index]);
// No more vector needed
if(rmax >= Rand)
{
for( iNumbQ2 = 1; iNumbQ2<length; iNumbQ2++ )
{
if(Rand <= pElD->TableCrossSec[index+iNumbQ2]) break;
}
}
/*
else // Build second part of the vector
{
if(!isIni)
{
R1 = pElD->R1;
R2 = pElD->R2;
Aeff = pElD->Aeff;
Pnucl = pElD->Pnucl;
hLabMomentum = std::sqrt(hLabMomentum2);
DefineHadronValues(Z);
Weight = pElD->CrossSecMaxQ2[NumbOnE];
}
G4int AWeight = pElD->AtomicWeight;
// Stop building when find out the node
for(iNumbQ2 = length; iNumbQ2<ONQ2; iNumbQ2++)
if(verboseLevel >0)
{
G4cout<<"1 plab T "<<plab<<" "<<T<<" sigTot B ReIm "
<<HadrTot<<" "<<HadrSlope<<" "<<HadrReIm<<G4endl;
G4cout<<" R1 R2 Aeff p "<<R1<<" "<<R2<<" "<<Aeff<<" "
<<Pnucl<<G4endl;
}
pElD->CrossSecMaxQ2[NumbOnE] = 1.0;
if(verboseLevel > 1)
G4cout<<" HadrNucleusQ2_2: NumbOnE= " << NumbOnE
<< " length= " << length
<< " Q2max= " << Q2max
<< " ekin= " << ekin <<G4endl;
pElD->TableCrossSec[index] = 0;
dQ2 = pElD->TableQ2[1]-pElD->TableQ2[0];
GetHeavyFq2(NumbN, LineFq2); // %%%%%%%%%%%%%%%%%%%%%%%%%
for(G4int ii=0; ii<ONQ2; ii++)
{
//if(verboseLevel > 2)
// G4cout<<" ii LineFq2 "<<ii<<" "<<LineFq2[ii]/LineFq2[ONQ2-1]
// <<" dF(q2) "<<LineFq2[ii]-LineFq2[ii-1]<<G4endl;
pElD->TableCrossSec[index+ii] = LineFq2[ii]/LineFq2[ONQ2-1];
}
pElD->dnkE[NumbOnE] = ONQ2;
length = ONQ2;
}
G4double* dNumbFQ2 = &(pElD->TableCrossSec[index]);
for( iNumbQ2 = 1; iNumbQ2<length; iNumbQ2++ )
{
Q2 = pElD->TableQ2[iNumbQ2];
if(Q2 < Q2max) Buf = GetLightFq2(Z, AWeight, Q2)/Weight;
else Buf = 1.0;
pElD->TableCrossSec[index+iNumbQ2] = Buf;
if(verboseLevel > 1)
G4cout<<" HadrNucleusQ2_2: NumbOnE= " << NumbOnE
<< " iNumbQ2= " << iNumbQ2 << " Q2= "
<<Q2 <<" Buf= " <<Buf<<" B*W "<<Buf*Weight<<G4endl;
if(Rand <= Buf)
{
pElD->dnkE[NumbOnE] = iNumbQ2+1;
break;
}
}
}
*/
if(Rand <= pElD->TableCrossSec[index+iNumbQ2]) break;
}
Q2 = GetQ2_2(iNumbQ2, dNumbQ2, dNumbFQ2, Rand);
if(tmax < Q2max) Q2 *= tmax/Q2max;
@@ -787,31 +685,32 @@ G4double G4ElasticHadrNucleusHE::GetQ2_2(G4int kk, G4double * Q,
<< D0 << G4endl;
if(std::abs(D0) < 0.00000001)
{
{
ranQ2 = X2 + (ranUni - F2)*(X3 - X2)/(F3 - F2);
}
}
else
{
{
G4double DA = X1*F2+X3*F1+X2*F3-X3*F2-X1*F3-X2*F1;
G4double DB = X2*F12+X1*F32+X3*F22-X2*F32-X3*F12-X1*F22;
G4double DC = X3*F2*F12+X2*F1*F32+X1*F3*F22
-X1*F2*F32-X2*F3*F12-X3*F1*F22;
ranQ2 = (DA*ranUni*ranUni + DB*ranUni + DC)/D0;
}
}
return ranQ2; // MeV^2
}
////////////////////////////////////////////////////////////////////////
//
//
G4double G4ElasticHadrNucleusHE::GetHeavyFq2(G4int Nucleus, G4double * LineF)
G4double G4ElasticHadrNucleusHE::GetHeavyFq2(G4int Nucleus, G4double* LineF)
{
G4int ii, jj, aSimp;
G4double Q2l, curQ2, ddQ2, curSum, curSec, totSum;
G4double curQ2, curSec;
G4double curSum = 0.0;
G4double totSum = 0.0;
ddQ2 = dQ2/20;
Q2l = 0;
totSum = 0;
G4double ddQ2 = dQ2/20;
G4double Q2l = 0;
LineF[0] = 0;
for(ii = 1; ii<ONQ2; ii++)
@@ -846,8 +745,9 @@ G4double G4ElasticHadrNucleusHE::GetHeavyFq2(G4int Nucleus, G4double * LineF)
<<curSec<<" totSum "<< totSum<<" DTot "
<<curSum<<G4endl;
}
return curSum;
return totSum;
}
////////////////////////////////////////////////////////////////////////
//
//
@@ -955,45 +855,7 @@ G4double G4ElasticHadrNucleusHE::GetLightFq2(G4int Z, G4int Nucleus,
if(std::fabs(N1*Prod1/Prod0) < prec) break;
} // i1
/*
for(G4int i1 = 1; i1<= Nucleus; i1++)
{
N1 *= UnucRho2*G4double(Nucleus-i1+1)/G4double(i1);
Prod1 = 0;
Tot0 = 0;
N2 = -1;
for(G4int i2 = 1; i2<=Nucleus; i2++)
{
N2 *= UnucRho2*G4double(Nucleus-i2+1)/G4double(i2);
Prod2 = 0;
N5 = -1.0/NN2;
for(G4int m2=0; m2<= i2; m2++)
{
Prod3 = 0;
exp2 = 1.0/(m2/R22B+(i2-m2)/R12B);
N5 *= (-NN2);
N4 = -1.0/NN2;
for(G4int m1=0; m1<=i1; m1++)
{
exp1 = 1.0/(m1/R22B+(i1-m1)/R12B);
dddd = 0.25*(exp1 + exp2);
N4 *= (-NN2);
Prod3 += N4*exp1*exp2*SetBinom[i1][m1]*(1-std::exp(-Q2*dddd))/dddd;
} // m1
Prod2 += Prod3*N5*SetBinom[i2][m2];
} // m2
Prod1 += Prod2*N2*std::cos(FiH*(i1-i2));
if (std::abs(Prod2*N2/Prod1)<prec) break;
} // i2
Prod0 += Prod1*N1;
if(std::abs(N1*Prod1/Prod0) < prec) break;
}
*/ // i1
Prod0 *= 0.25*pi/MbToGeV2; // This is in mb
if(verboseLevel>1)
@@ -1002,13 +864,13 @@ G4double G4ElasticHadrNucleusHE::GetLightFq2(G4int Z, G4int Nucleus,
return Prod0;
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
G4double G4ElasticHadrNucleusHE::
G4double G4ElasticHadrNucleusHE::
HadrNucDifferCrSec(G4int Nucleus, G4double aQ2)
{
// ------ All external kinematical variables are in MeV -------
// ------ but internal in GeV !!!! ------
G4double theQ2 = aQ2; ///GeV/GeV;
G4double theQ2 = aQ2; ///GeV/GeV;
// Scattering of proton
if(Nucleus == 1)
@@ -1149,41 +1011,41 @@ G4double G4ElasticHadrNucleusHE::GetLightFq2(G4int Z, G4int Nucleus,
C2/exp2p*std::exp(-theQ2/4/exp2p)+
C3/exp3p*std::exp(-theQ2/4/exp3p));
DmedTot = DmedTot + N2p*BinCoeff*
(C1/exp1p-C2/exp2p+C3/exp3p);
DmedTot = DmedTot + N2p*BinCoeff*
(C1/exp1p-C2/exp2p+C3/exp3p);
N2p = -N2p*R23dR13;
} // l
N2p = -N2p*R23dR13;
} // l
Din1 = Din1+Din2*N1p/*Mnoj[i]*//(i+2)/(i+1)*std::cos(FiH*i);
DTot1 = DTot1+DmedTot*N1p/*Mnoj[i]*//(i+2)/(i+1)*std::cos(FiH*i);
if(std::fabs(Din2*N1p/Din1) < 0.000001) break;
if(std::fabs(Din2*N1p/Din1) < 0.000001) break;
} // i
Din1 = -Din1*Nucleus*(Nucleus-1)
Din1 = -Din1*Nucleus*(Nucleus-1)
/2/pi/Normp/2/pi/Normp*16*pi*pi;
DTot1 = DTot1*Nucleus*(Nucleus-1)
DTot1 = DTot1*Nucleus*(Nucleus-1)
/2/pi/Normp/2/pi/Normp*16*pi*pi;
DTot1 *= 5; // $$$$$$$$$$$$$$$$$$$$$$$$
DTot1 *= 5; // $$$$$$$$$$$$$$$$$$$$$$$$
// Din1 *= 0.2; // %%%%%%%%%%%%%%%%%%%%%%% proton
// Din1 *= 0.05; // %%%%%%%%%%%%%%%%%%%%%%% pi+
// ---------------- dSigma/d|-t|, mb/(GeV/c)^-2 -----------------
G4double DiffCrSec2 = (ReElasticAmpl0*ReElasticAmpl0+
G4double DiffCrSec2 = (ReElasticAmpl0*ReElasticAmpl0+
(ImElasticAmpl0+Din1)*
(ImElasticAmpl0+Din1))*2/4/pi;
Tot1 = Tot1-DTot1;
Tot1 = Tot1-DTot1;
// Tott1 = Tot1*1.0;
Dtot11 = DTot1;
aAIm = ImElasticAmpl0;
aDIm = Din1;
Dtot11 = DTot1;
aAIm = ImElasticAmpl0;
aDIm = Din1;
return DiffCrSec2*1.0; // dSig/d|-t|, mb/(GeV/c)^-2
} // function
return DiffCrSec2*1.0; // dSig/d|-t|, mb/(GeV/c)^-2
} // function
// ##############################################
////////////////////////////////////////////////////////////////
@@ -1199,7 +1061,7 @@ void G4ElasticHadrNucleusHE::DefineHadronValues(G4int Z)
G4double Ecm = 0.5*(sHadr-hMass2+protonM2)/sqrS;
MomentumCM = std::sqrt(Ecm*Ecm-protonM2);
if(verboseLevel>2)
if(verboseLevel>2)
G4cout << "GetHadrVall.: Z= " << Z << " iHadr= " << iHadron
<< " E(GeV)= " << HadrEnergy << " sqrS= " << sqrS
<< " plab= " << hLabMomentum
@@ -1216,52 +1078,52 @@ void G4ElasticHadrNucleusHE::DefineHadronValues(G4int Z)
case 6:
if(hLabMomentum > 10)
TotP = TotN = 7.5*logE - 40.12525 + 103*std::pow(sHadr,-0.165); // mb
TotP = TotN = 7.5*logE - 40.12525 + 103*std::pow(sHadr,-0.165); // mb
else
{
{
// ================== neutron ================
//// if(iHadrCode == 2112)
if( hLabMomentum > 1.4 )
TotN = 33.3+15.2*(hLabMomentum2-1.35)/
(std::pow(hLabMomentum,2.37)+0.95);
if( hLabMomentum > 1.4 )
TotN = 33.3+15.2*(hLabMomentum2-1.35)/
(std::pow(hLabMomentum,2.37)+0.95);
else if(hLabMomentum > 0.8)
{
G4double A0 = logE + 0.0513;
TotN = 33.0 + 25.5*A0*A0;
}
else
{
G4double A0 = logE - 0.2634; // log(1.3)
TotN = 33.0 + 30.*A0*A0*A0*A0;
}
else if(hLabMomentum > 0.8)
{
G4double A0 = logE + 0.0513;
TotN = 33.0 + 25.5*A0*A0;
}
else
{
G4double A0 = logE - 0.2634; // log(1.3)
TotN = 33.0 + 30.*A0*A0*A0*A0;
}
// ================= proton ===============
// else if(iHadrCode == 2212)
{
if(hLabMomentum >= 1.05)
{
if(hLabMomentum >= 1.05)
{
TotP = 39.0+75.*(hLabMomentum-1.2)/
(hLabMomentum2*hLabMomentum+0.15);
TotP = 39.0+75.*(hLabMomentum-1.2)/
(hLabMomentum2*hLabMomentum+0.15);
}
else if(hLabMomentum >= 0.7)
else if(hLabMomentum >= 0.7)
{
G4double A0 = logE + 0.3147;
TotP = 23.0 + 40.*A0*A0;
}
else
else
{
TotP = 23.+50.*std::pow(std::log(0.73/hLabMomentum),3.5);
TotP = 23.+50.*std::pow(std::log(0.73/hLabMomentum),3.5);
}
}
}
}
}
// HadrTot = 0.5*(82*TotP+126*TotN)/104; // $$$$$$$$$$$$$$$$$$
HadrTot = 0.5*(TotP+TotN);
HadrTot = 0.5*(TotP+TotN);
// ...................................................
// Proton slope
if(hLabMomentum >= 2.) HadrSlope = 5.44 + 0.88*logS;
@@ -1276,7 +1138,7 @@ void G4ElasticHadrNucleusHE::DefineHadronValues(G4int Z)
else if(hLabMomentum >= 0.6)
HadrReIm = -75.5*(std::pow(hLabMomentum,0.25)-0.95)/
(std::pow(3*hLabMomentum,2.2)+1);
(std::pow(3*hLabMomentum,2.2)+1);
else
HadrReIm = 15.5*hLabMomentum/(27*hLabMomentum2*hLabMomentum+2);
@@ -1449,10 +1311,6 @@ void G4ElasticHadrNucleusHE::DefineHadronValues(G4int Z)
<< " HadrReIm= " << HadrReIm << " DDSect2= " << DDSect2
<< " DDSect3= " << DDSect3 << G4endl;
// HadrTot1 = HadrTot;
// HadrSlope1 = HadrSlope;
// HadrReIm1 = HadrReIm;
if(Z != 1) return;
// Scattering of protons
@@ -1511,10 +1369,9 @@ void G4ElasticHadrNucleusHE::DefineHadronValues(G4int Z)
static const G4double B1KM[2]={1.6,1.6};
switch(iHadron)
{
{
case 0 :
if(hLabMomentum <BoundaryP[0])
InterpolateHN(6,EnP0,C0P0,C1P0,B0P0,B1P0);
@@ -1522,7 +1379,6 @@ void G4ElasticHadrNucleusHE::DefineHadronValues(G4int Z)
break;
case 6 :
// case 7 :
if(hLabMomentum < BoundaryP[1])
InterpolateHN(5,EnN,C0N,C1N,B0N,B1N);
@@ -1568,95 +1424,83 @@ void G4ElasticHadrNucleusHE::DefineHadronValues(G4int Z)
if(hLabMomentum < 1) Coeff2 = 0.01;
else Coeff2 = 0.01/hLabMomentum2/hLabMomentum;
break;
}
}
if(verboseLevel > 2)
G4cout<<" HadrVal : Plasb "<<hLabMomentum
<<" iHadron "<<iHadron<<" HadrTot "<<HadrTot<<G4endl;
if(verboseLevel > 2)
G4cout<<" HadrVal : Plasb "<<hLabMomentum
<<" iHadron "<<iHadron<<" HadrTot "<<HadrTot<<G4endl;
}
// =====================================================
void G4ElasticHadrNucleusHE::
void G4ElasticHadrNucleusHE::
GetKinematics(const G4ParticleDefinition * aHadron,
G4double MomentumH)
{
{
if (verboseLevel>1)
G4cout<<"1 GetKin.: HadronName MomentumH "
<<aHadron->GetParticleName()<<" "<<MomentumH<<G4endl;
G4cout<<"1 GetKin.: HadronName MomentumH "
<<aHadron->GetParticleName()<<" "<<MomentumH<<G4endl;
DefineHadronValues(1);
DefineHadronValues(1);
G4double Sh = 2.0*protonM*HadrEnergy+protonM2+hMass2; // GeV
G4double Sh = 2.0*protonM*HadrEnergy+protonM2+hMass2; // GeV
ConstU = 2*protonM2+2*hMass2-Sh;
ConstU = 2*protonM2+2*hMass2-Sh;
G4double MaxT = 4*MomentumCM*MomentumCM;
G4double MaxT = 4*MomentumCM*MomentumCM;
BoundaryTL[0] = MaxT; //2.0;
BoundaryTL[1] = MaxT;
BoundaryTL[3] = MaxT;
BoundaryTL[4] = MaxT;
BoundaryTL[5] = MaxT;
BoundaryTL[0] = MaxT; //2.0;
BoundaryTL[1] = MaxT;
BoundaryTL[3] = MaxT;
BoundaryTL[4] = MaxT;
BoundaryTL[5] = MaxT;
G4int NumberH=0;
G4int NumberH=0;
while(iHadrCode!=HadronCode[NumberH]) NumberH++;
while(iHadrCode!=HadronCode[NumberH]) NumberH++;
NumberH = HadronType1[NumberH];
NumberH = HadronType1[NumberH];
if(MomentumH<BoundaryP[NumberH]) MaxTR = BoundaryTL[NumberH];
else MaxTR = BoundaryTG[NumberH];
if(MomentumH<BoundaryP[NumberH]) MaxTR = BoundaryTL[NumberH];
else MaxTR = BoundaryTG[NumberH];
if (verboseLevel>1)
G4cout<<"3 GetKin. : NumberH "<<NumberH
<<" Bound.P[NumberH] "<<BoundaryP[NumberH]
<<" Bound.TL[NumberH] "<<BoundaryTL[NumberH]
<<" Bound.TG[NumberH] "<<BoundaryTG[NumberH]
<<" MaxT MaxTR "<<MaxT<<" "<<MaxTR<<G4endl;
G4cout<<"3 GetKin. : NumberH "<<NumberH
<<" Bound.P[NumberH] "<<BoundaryP[NumberH]
<<" Bound.TL[NumberH] "<<BoundaryTL[NumberH]
<<" Bound.TG[NumberH] "<<BoundaryTG[NumberH]
<<" MaxT MaxTR "<<MaxT<<" "<<MaxTR<<G4endl;
// GetParametersHP(aHadron, MomentumH);
}
}
// ============================================================
G4double G4ElasticHadrNucleusHE::GetFt(G4double Q2)
{
G4float Fdistr=0, SqrQ2 = std::sqrt(Q2);
G4double G4ElasticHadrNucleusHE::GetFt(G4double Q2)
{
G4double Fdistr=0;
G4double SqrQ2 = std::sqrt(Q2);
Fdistr = (1-Coeff1-Coeff0) //-0.0*Coeff2*std::exp(ConstU))
/HadrSlope*(1-std::exp(-HadrSlope*Q2))
+Coeff0*(1-std::exp(-Slope0*Q2))
+Coeff2/Slope2*std::exp(Slope2*ConstU)*
(std::exp(Slope2*Q2)-1)
+2*Coeff1/Slope1*(1/Slope1-(1/Slope1+SqrQ2)*
std::exp(-Slope1*SqrQ2))
;
Fdistr = (1-Coeff1-Coeff0) //-0.0*Coeff2*std::exp(ConstU))
/HadrSlope*(1-std::exp(-HadrSlope*Q2))
+ Coeff0*(1-std::exp(-Slope0*Q2))
+ Coeff2/Slope2*std::exp(Slope2*ConstU)*(std::exp(Slope2*Q2)-1)
+ 2*Coeff1/Slope1*(1/Slope1-(1/Slope1+SqrQ2)*std::exp(-Slope1*SqrQ2));
if (verboseLevel>1)
G4cout<<"Old: Coeff0 Coeff1 Coeff2 "<<Coeff0<<" "
<<Coeff1<<" "<<Coeff2<<" Slope Slope0 Slope1 Slope2 "
<<HadrSlope<<" "<<Slope0<<" "<<Slope1<<" "<<Slope2
<<" Fdistr "<<Fdistr<<G4endl;
return Fdistr;
}
return Fdistr;
}
// +++++++++++++++++++++++++++++++++++++++
G4double G4ElasticHadrNucleusHE::
GetDistrFun(G4double Q2)
{
return GetFt(Q2)/FmaxT;
}
// +++++++++++++++++++++++++++++++++++++++
G4double G4ElasticHadrNucleusHE::
GetQ2(G4double Ran)
{
G4double DDD0=MaxTR*0.5, DDD1=0.0, DDD2=MaxTR, delta;
G4double Q2=0;
G4double G4ElasticHadrNucleusHE::GetQ2(G4double Ran)
{
G4double DDD0=MaxTR*0.5, DDD1=0.0, DDD2=MaxTR, delta;
G4double Q2=0;
FmaxT = GetFt(MaxTR);
delta = GetDistrFun(DDD0)-Ran;
FmaxT = GetFt(MaxTR);
delta = GetDistrFun(DDD0)-Ran;
while(std::fabs(delta) > 0.0001)
while(std::fabs(delta) > 0.0001)
{
if(delta>0)
{
@@ -1668,34 +1512,33 @@ if(verboseLevel > 2)
DDD1 = DDD0;
DDD0 = (DDD0+DDD2)*0.5;
}
delta = GetDistrFun(DDD0)-Ran;
delta = GetDistrFun(DDD0)-Ran;
}
Q2 = DDD0;
Q2 = DDD0;
return Q2;
}
return Q2;
}
// ++++++++++++++++++++++++++++++++++++++++++
G4double G4ElasticHadrNucleusHE::
G4double G4ElasticHadrNucleusHE::
HadronProtonQ2(const G4ParticleDefinition * p,
G4double inLabMom)
{
{
hMass = p->GetPDGMass()/GeV;
hMass2 = hMass*hMass;
hLabMomentum = inLabMom;
hLabMomentum2 = hLabMomentum*hLabMomentum;
HadrEnergy = sqrt(hLabMomentum2+hMass2);
hMass = p->GetPDGMass()/GeV;
hMass2 = hMass*hMass;
hLabMomentum = inLabMom;
hLabMomentum2 = hLabMomentum*hLabMomentum;
HadrEnergy = sqrt(hLabMomentum2+hMass2);
G4double Rand = G4UniformRand();
G4double Rand = G4UniformRand(), Q2=0;
GetKinematics(p, inLabMom);
GetKinematics(p, inLabMom);
G4double Q2 = GetQ2(Rand);
Q2 = GetQ2(Rand);
return Q2;
}
return Q2;
}
// ===========================================
source/processes/hadronic/models/coherent_elastic/src/G4HadronElastic.cc
+86 -51
View File
@@ -23,8 +23,8 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4HadronElastic.cc,v 1.55 2007/12/10 09:56:01 vnivanch Exp $
// GEANT4 tag $Name: geant4-09-01 $
// $Id: G4HadronElastic.cc,v 1.61 2008/08/05 07:37:39 vnivanch Exp $
// GEANT4 tag $Name: geant4-09-02 $
//
//
// Physics model class G4HadronElastic (derived from G4LElastic)
@@ -86,8 +86,8 @@ G4HadronElastic::G4HadronElastic(G4ElasticHadrNucleusHE* HModel)
verboseLevel= 0;
lowEnergyRecoilLimit = 100.*keV;
lowEnergyLimitQ = 0.0*GeV;
lowEnergyLimitHE = DBL_MAX;
lowestEnergyLimit= 0.0*keV;
lowEnergyLimitHE = 1.0*GeV;
lowestEnergyLimit= 1.e-6*eV;
plabLowLimit = 20.0*MeV;
qCManager = G4QElasticCrossSection::GetPointer();
@@ -99,11 +99,24 @@ G4HadronElastic::G4HadronElastic(G4ElasticHadrNucleusHE* HModel)
theAlpha = G4Alpha::Alpha();
thePionPlus = G4PionPlus::PionPlus();
thePionMinus= G4PionMinus::PionMinus();
nnans = 0;
npos = 0;
nneg = 0;
neneg = 0;
}
G4HadronElastic::~G4HadronElastic()
{
delete hElastic;
if( (nnans + npos + nneg + neneg) > 0 ) {
G4cout << "### G4HadronElastic destructor Warnings: ";
if(nnans > 0) G4cout << "### N(nans) = " << nnans;
if(npos > 0) G4cout << "### N(cost > 1)= " << npos;
if(nneg > 0) G4cout << "### N(cost <-1)= " << nneg;
if(neneg > 0) G4cout << "### N(E < 0)= " << neneg;
G4cout << "###" << G4endl;
}
}
G4VQCrossSection* G4HadronElastic::GetCS()
@@ -133,12 +146,12 @@ G4HadFinalState* G4HadronElastic::ApplyYourself(
G4double zTarget = targetNucleus.GetZ();
G4double plab = aParticle->GetTotalMomentum();
if (verboseLevel >1)
if (verboseLevel >1) {
G4cout << "G4HadronElastic::DoIt: Incident particle plab="
<< plab/GeV << " GeV/c "
<< " ekin(MeV) = " << ekin/MeV << " "
<< aParticle->GetDefinition()->GetParticleName() << G4endl;
}
// Scattered particle referred to axis of incident particle
const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
G4double m1 = theParticle->GetPDGMass();
@@ -147,12 +160,12 @@ G4HadFinalState* G4HadronElastic::ApplyYourself(
G4int A = static_cast<G4int>(aTarget+0.5);
G4int N = A - Z;
G4int projPDG = theParticle->GetPDGEncoding();
if (verboseLevel>1)
if (verboseLevel>1) {
G4cout << "G4HadronElastic for " << theParticle->GetParticleName()
<< " PDGcode= " << projPDG << " on nucleus Z= " << Z
<< " A= " << A << " N= " << N
<< G4endl;
}
G4ParticleDefinition * theDef = 0;
if(Z == 1 && A == 1) theDef = theProton;
@@ -179,27 +192,27 @@ G4HadFinalState* G4HadronElastic::ApplyYourself(
G4ElasticGenerator gtype = fLElastic;
// Q-elastic for p,n scattering on H and He
if (theParticle == theProton || theParticle == theNeutron)
if (theParticle == theProton || theParticle == theNeutron) {
// && Z <= 2 && ekin >= lowEnergyLimitQ)
gtype = fQElastic;
else {
} else {
// S-wave for very low energy
if(plab < plabLowLimit) gtype = fSWave;
// HE-elastic for energetic projectile mesons
// else if(ekin >= lowEnergyLimitHE && theParticle->GetBaryonNumber() == 0)
else if(ekin >= lowEnergyLimitHE && (theParticle == thePionPlus || theParticle == thePionMinus))
gtype = fHElastic;
else if(ekin >= lowEnergyLimitHE && theParticle->GetBaryonNumber() == 0)
{ gtype = fHElastic; }
}
//
// Sample t
//
if(gtype == fQElastic) {
if (verboseLevel >1)
if (verboseLevel >1) {
G4cout << "G4HadronElastic: Z= " << Z << " N= "
<< N << " pdg= " << projPDG
<< " mom(GeV)= " << plab/GeV << " " << qCManager << G4endl;
}
if(Z == 1 && N == 2) N = 1;
else if(Z == 2 && N == 1) N = 2;
G4double cs = qCManager->GetCrossSection(false,plab,Z,N,projPDG);
@@ -211,7 +224,8 @@ G4HadFinalState* G4HadronElastic::ApplyYourself(
}
if(gtype == fLElastic) {
t = GeV*GeV*SampleT(ptot,m1,m2,aTarget);
G4double g2 = GeV*GeV;
t = g2*SampleT(tmax/g2,m1,m2,aTarget);
}
// use mean atomic number
@@ -233,36 +247,48 @@ G4HadFinalState* G4HadronElastic::ApplyYourself(
<< G4endl;
}
t = 0.0;
nnans++;
}
if(gtype == fSWave) t = G4UniformRand()*tmax;
if(verboseLevel>1)
if(verboseLevel>1) {
G4cout <<"type= " << gtype <<" t= " << t << " tmax= " << tmax
<< " ptot= " << ptot << G4endl;
}
// Sampling in CM system
G4double phi = G4UniformRand()*twopi;
G4double cost = 1. - 2.0*t/tmax;
G4double sint;
if( cost >= 1.0 )
{
// problem in sampling
if(cost >= 1.0) {
cost = 1.0;
sint = 0.0;
}
else if( cost <= -1.0)
{
cost = -1.0;
sint = 0.0;
}
else
{
npos++;
} else if(cost < -1 ) {
/*
G4cout << "G4HadronElastic:WARNING: Z= " << Z << " N= "
<< N << " " << aParticle->GetDefinition()->GetParticleName()
<< " mom(GeV)= " << plab/GeV
<< " the model type " << gtype;
if(gtype == fQElastic) G4cout << " CHIPS ";
else if(gtype == fLElastic) G4cout << " LElastic ";
else if(gtype == fHElastic) G4cout << " HElastic ";
G4cout << " cost= " << cost
<< G4endl;
*/
cost = 1.0;
sint = 0.0;
nneg++;
// normal situation
} else {
sint = std::sqrt((1.0-cost)*(1.0+cost));
}
if (verboseLevel>1)
if (verboseLevel>1) {
G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
}
G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
v1 *= ptot;
G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
@@ -270,32 +296,37 @@ G4HadFinalState* G4HadronElastic::ApplyYourself(
nlv1.boost(bst);
G4double eFinal = nlv1.e() - m1;
if (verboseLevel > 1)
if (verboseLevel > 1) {
G4cout << "Scattered: "
<< nlv1<<" m= " << m1 << " ekin(MeV)= " << eFinal
<< " Proj: 4-mom " << lv1
<<G4endl;
if(eFinal < 0.0) {
G4cout << "G4HadronElastic WARNING ekin= " << eFinal
<< " after scattering of "
<< aParticle->GetDefinition()->GetParticleName()
<< " p(GeV/c)= " << plab
<< " on " << theDef->GetParticleName()
<< G4endl;
eFinal = 0.0;
nlv1.setE(m1);
}
if(eFinal <= lowestEnergyLimit) {
if(eFinal < 0.0 && verboseLevel > 0) {
neneg++;
G4cout << "G4HadronElastic WARNING ekin= " << eFinal
<< " after scattering of "
<< aParticle->GetDefinition()->GetParticleName()
<< " p(GeV/c)= " << plab
<< " on " << theDef->GetParticleName()
<< G4endl;
}
theParticleChange.SetEnergyChange(0.0);
nlv1 = G4LorentzVector(0.0,0.0,0.0,m1);
} else {
theParticleChange.SetMomentumChange(nlv1.vect().unit());
theParticleChange.SetEnergyChange(eFinal);
}
theParticleChange.SetMomentumChange(nlv1.vect().unit());
theParticleChange.SetEnergyChange(eFinal);
G4LorentzVector nlv0 = lv - nlv1;
G4double erec = nlv0.e() - m2;
if (verboseLevel > 1)
if (verboseLevel > 1) {
G4cout << "Recoil: "
<< nlv0<<" m= " << m2 << " ekin(MeV)= " << erec
<<G4endl;
}
if(erec > lowEnergyRecoilLimit) {
G4DynamicParticle * aSec = new G4DynamicParticle(theDef, nlv0);
theParticleChange.AddSecondary(aSec);
@@ -308,14 +339,14 @@ G4HadFinalState* G4HadronElastic::ApplyYourself(
}
G4double
G4HadronElastic::SampleT(G4double, G4double, G4double, G4double atno2)
G4HadronElastic::SampleT(G4double tmax, G4double, G4double, G4double atno2)
{
// G4cout << "Entering elastic scattering 2"<<G4endl;
// Compute the direction of elastic scattering.
// It is planned to replace this code with a method based on
// parameterized functions and a Monte Carlo method to invert the CDF.
G4double ran = G4UniformRand();
// G4double ran = G4UniformRand();
G4double aa, bb, cc, dd, rr;
if (atno2 <= 62.) {
aa = std::pow(atno2, 1.63);
@@ -330,14 +361,18 @@ G4HadronElastic::SampleT(G4double, G4double, G4double, G4double atno2)
}
aa = aa/bb;
cc = cc/dd;
G4double ran, t1, t2;
do {
ran = G4UniformRand();
t1 = -std::log(ran)/bb;
t2 = -std::log(ran)/dd;
} while(t1 > tmax || t2 > tmax);
rr = (aa + cc)*ran;
if (verboseLevel > 1) {
G4cout << "DoIt: aa,bb,cc,dd,rr" << G4endl;
G4cout << aa << " " << bb << " " << cc << " " << dd << " " << rr << G4endl;
}
G4double t1 = -std::log(ran)/bb;
G4double t2 = -std::log(ran)/dd;
if (verboseLevel > 1) {
G4cout << "t1,Fctcos " << t1 << " " << Fctcos(t1, aa, bb, cc, dd, rr) << G4endl;
G4cout << "t2,Fctcos " << t2 << " " << Fctcos(t2, aa, bb, cc, dd, rr) << G4endl;
}
@@ -348,7 +383,7 @@ G4HadronElastic::SampleT(G4double, G4double, G4double, G4double atno2)
ier1 = Rtmi(&t, t1, t2, eps, ind1,
aa, bb, cc, dd, rr);
if (verboseLevel > 1) {
G4cout << "From Rtmi, ier1=" << ier1 << G4endl;
G4cout << "From Rtmi, ier1=" << ier1 << " t= " << t << G4endl;
G4cout << "t, Fctcos " << t << " " << Fctcos(t, aa, bb, cc, dd, rr) << G4endl;
}
if (ier1 != 0) t = 0.25*(3.*t1 + t2);
source/processes/hadronic/models/coherent_elastic/src/G4UHadronElasticProcess.cc
+13 -9
View File
@@ -23,8 +23,8 @@
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: G4UHadronElasticProcess.cc,v 1.35 2007/11/15 15:53:25 vnivanch Exp $
// GEANT4 tag $Name: geant4-09-01 $
// $Id: G4UHadronElasticProcess.cc,v 1.39 2008/10/22 08:16:40 vnivanch Exp $
// GEANT4 tag $Name: geant4-09-02 $
//
// Geant4 Hadron Elastic Scattering Process -- header file
//
@@ -57,6 +57,7 @@
G4UHadronElasticProcess::G4UHadronElasticProcess(const G4String& pName, G4double)
: G4HadronicProcess(pName), lowestEnergy(0.0), first(true)
{
SetProcessSubType(fHadronElastic);
AddDataSet(new G4HadronElasticDataSet);
theProton = G4Proton::Proton();
theNeutron = G4Neutron::Neutron();
@@ -88,9 +89,10 @@ BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
// defined lowest threshold for the projectile
if(theParticle->GetPDGCharge() != 0.0) lowestEnergy = eV;
if(verboseLevel>1 ||
(verboseLevel==1 && theParticle == theNeutron)) {
G4cout << G4endl;
// if(verboseLevel>1 ||
// (verboseLevel==1 && theParticle == theNeutron)) {
if(verboseLevel>1 && theParticle == theNeutron) {
// G4cout << G4endl;
G4cout << "G4UHadronElasticProcess for "
<< theParticle->GetParticleName()
<< " PDGcode= " << pPDG
@@ -99,7 +101,8 @@ BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
<< G4endl;
}
}
store->BuildPhysicsTable(aParticleType);
G4HadronicProcess::BuildPhysicsTable(aParticleType);
//store->BuildPhysicsTable(aParticleType);
}
G4double G4UHadronElasticProcess::GetMeanFreePath(const G4Track& track,
@@ -168,7 +171,7 @@ G4double G4UHadronElasticProcess::GetMicroscopicCrossSection(
x = 0.0;
if(ni == 0) {
G4int N = G4int(elm->GetN()+0.5) - iz;
x = qCManager->GetCrossSection(true,momentum,iz,N,pPDG);
x = qCManager->GetCrossSection(false,momentum,iz,N,pPDG);
xsecH[0] = x;
#ifdef G4VERBOSE
if(verboseLevel>1)
@@ -177,7 +180,6 @@ G4double G4UHadronElasticProcess::GetMicroscopicCrossSection(
<< " mom(GeV)= " << momentum/GeV
<< " " << qCManager << G4endl;
#endif
} else {
G4double* ab = elm->GetRelativeAbundanceVector();
for(G4int j=0; j<ni; j++) {
@@ -187,12 +189,14 @@ G4double G4UHadronElasticProcess::GetMicroscopicCrossSection(
} else {
N = 2;
}
#ifdef G4VERBOSE
if(verboseLevel>1)
G4cout << "G4UHadronElasticProcess compute CHIPS CS for Z= " << iz
<< " N= " << N << " pdg= " << pPDG
<< " mom(GeV)= " << momentum/GeV
<< " " << qCManager << G4endl;
G4double y = ab[j]*qCManager->GetCrossSection(true,momentum,iz,N,pPDG);
#endif
G4double y = ab[j]*qCManager->GetCrossSection(false,momentum,iz,N,pPDG);
x += y;
xsecH[j] = x;
}