61 : myA(0), myZ(0), myL(0), theNucleons(250), currentNucleon(-1), theDensity(0),
62 nucleondistance(0.8*
fermi),excitationEnergy(0.),
63 places(250), momentum(250), fermiM(250), testSums(250)
72#if defined(NON_INTEGER_A_Z)
123 for (
G4int aNucleon=0; aNucleon <
myA; aNucleon++)
193 for (
int i=0; i<
myA; i++)
195 if (
theNucleons[i].GetPosition().mag2() > maxradius2 )
231 G4double factor=(1-std::sqrt(1-beta2))/beta2;
235 factor * (theBeta*
theNucleons[i].GetPosition()) * theBeta;
243 if (theBoost.
e() !=0 ) {
282 G4int protons=0, nucleons=0, lambdas=0;
285 while ( nucleons <
myA ) {
287 if ( rnd < probProton ) {
288 if ( protons <
myZ ) {
292 }
else if ( rnd < probProton + probLambda ) {
293 if ( lambdas <
myL ) {
298 if ( (nucleons - protons - lambdas) < (
myA -
myZ -
myL) ) {
322 while ( (i <
myA) && (--interationsLeft>0))
329 G4RandFlat::shootArray(jr,prand);
334 aPos.
set((2*arand[jx]-1.), (2*arand[jy]-1.), (2*arand[--jr]-1.));
335 }
while (aPos.
mag2() > 1. );
341 std::vector<G4ThreeVector>::iterator iplace;
342 for( iplace=
places.begin(); iplace!=
places.end() && freeplace;++iplace)
344 delta = *iplace - aPos;
345 freeplace= delta.
mag2() > nd2;
355 G4double eFermi= std::sqrt(
sqr(pFermi) +
sqr(nucMass) ) - nucMass;
366 if (interationsLeft<=0) {
368 "Problem to place nucleons");
385 G4int loopCounterLeft = 10000;
386 for(
G4int ii=1; ii<4; ii++)
394 for(
G4int jj=0; jj < ii; jj++)
396 if( (
theNucleons[ii].GetPosition() -
theNucleons[jj].GetPosition()).mag2() <= nd2 ) {Continue =
true;
break;}
398 }
while( Continue && --loopCounterLeft > 0 );
400 if ( loopCounterLeft <= 0 ) {
402 "Unable to find a good position for the first alpha cluster");
404 loopCounterLeft = 10000;
405 for(
G4int ii=4; ii<8; ii++)
413 for(
G4int jj=0; jj < ii; jj++)
415 if( (
theNucleons[ii].GetPosition() -
theNucleons[jj].GetPosition()).mag2() <= nd2 ) {Continue =
true;
break;}
417 }
while( Continue && --loopCounterLeft > 0 );
419 if ( loopCounterLeft <= 0 ) {
421 "Unable to find a good position for the second alpha cluster");
423 loopCounterLeft = 10000;
424 for(
G4int ii=8; ii<12; ii++)
432 for(
G4int jj=0; jj < ii; jj++)
434 if( (
theNucleons[ii].GetPosition() -
theNucleons[jj].GetPosition()).mag2() <= nd2 ) {Continue =
true;
break;}
436 }
while( Continue && --loopCounterLeft > 0 );
438 if ( loopCounterLeft <= 0 ) {
440 "Unable to find a good position for the third alpha cluster");
466 for (
G4int ntry=0; ntry<1 ; ntry ++ )
468 for (i=0; i <
myA; i++ )
477 if ( eMax >
theNucleons[i].GetDefinition()->GetPDGMass() )
480 fermiM[i] = std::sqrt(pmax2);
481 while ( mom.
mag2() > pmax2 )
491 ed <<
"proton with eMax=" << eMax <<
G4endl;
492 G4Exception(
"G4Fancy3DNucleus::ChooseFermiMomenta(): difficulty finding proton momentum, set it to (0,0,0)",
511 for ( i=0; i<
myA ; i++ )
533 if ( sum.
mag() <= PFermi )
545 for (
G4int aNucleon=0; aNucleon <
myA-1; aNucleon++) {
546 delta = 2.*((
momentum[aNucleon]*testDir)*testDir);
548 testSums[aNucleon].Fill(delta, delta.
mag(), aNucleon);
555 while ( (sum-
testSums[--index].Vector).mag()>PFermi && index>0)
558 if ( sum.
mag() > (sum-
testSums[index].Vector).mag() ) {
564 if ( (sum-
testSums[index].Vector).mag() <= PFermi )
568 for (
G4int aNucleon=0; aNucleon<=index; aNucleon++)
573 && std::abs(
momentum[
myA-1].mag() - pTry ) < pBest )
581 G4String text =
"G4Fancy3DNucleus.cc: Logic error in ReduceSum()";
595 if (
fermiM[++swapit] > PFermi )
break;
597 if (swapit ==
myA-1 )
return false;
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
bool G4Fancy3DNucleusHelperForSortInZ(const G4Nucleon &nuc1, const G4Nucleon &nuc2)
CLHEP::HepLorentzVector G4LorentzVector
static constexpr double fermi
static constexpr double GeV
static constexpr double MeV
static constexpr double pi
CLHEP::Hep3Vector G4ThreeVector
void set(double x, double y, double z)
HepLorentzRotation & rotateY(double delta)
HepLorentzRotation & rotateZ(double delta)
const std::vector< G4Nucleon > & GetNucleons()
std::vector< G4double > fermiM
std::vector< G4Nucleon > theNucleons
G4double CoulombBarrier()
std::vector< G4ThreeVector > momentum
G4Nucleon * GetNextNucleon()
std::vector< G4Fancy3DNucleusHelper > testSums
void ChooseFermiMomenta()
const G4VNuclearDensity * GetNuclearDensity() const
void DoLorentzContraction(const G4LorentzVector &theBoost)
void Init(G4int theA, G4int theZ, G4int numberOfLambdas=0)
G4double GetNuclearRadius()
void DoTranslation(const G4ThreeVector &theShift)
G4double GetOuterRadius()
G4VNuclearDensity * theDensity
void DoLorentzBoost(const G4LorentzVector &theBoost)
G4double excitationEnergy
std::vector< G4ThreeVector > places
G4ThreeVector GetMomentum(G4double density, G4double maxMomentum=-1.)
G4double GetFermiMomentum(G4double density)
void Init(G4int anA, G4int aZ)
static G4double GetNuclearMass(G4int A, G4int Z, G4int L)
static G4Lambda * Lambda()
static G4Neutron * Neutron()
static G4double GetBindingEnergy(const G4int A, const G4int Z)
const G4ThreeVector & GetPosition() const
G4double GetPDGMass() const
static G4Pow * GetInstance()
G4double Z13(G4int Z) const
static G4Proton * Proton()
G4double GetDensity(const G4ThreeVector &aPosition) const
virtual G4double GetRelativeDensity(const G4ThreeVector &aPosition) const =0
virtual G4double GetRadius(const G4double maxRelativeDenisty) const =0
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double energy(const ThreeVector &p, const G4double m)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments