60 : myA(0), myZ(0), theNucleons(250), currentNucleon(-1), theDensity(0),
61 nucleondistance(0.8*
fermi),excitationEnergy(0.),
62 places(250), momentum(250), fermiM(250), testSums(250)
71 #if defined(NON_INTEGER_A_Z)
118 for (
G4int aNucleon=0; aNucleon <
myA; aNucleon++)
188 for (
int i=0; i<
myA; i++)
190 if (
theNucleons[i].GetPosition().mag2() > maxradius2 )
225 G4double factor=(1-std::sqrt(1-beta2))/beta2;
229 factor * (theBeta*
theNucleons[i].GetPosition()) * theBeta;
237 if (theBoost.
e() !=0 ) {
276 G4int protons=0,nucleons=0;
278 while (nucleons <
myA )
285 else if ( (nucleons-protons) < (
myA-
myZ) )
289 else G4cout <<
"G4Fancy3DNucleus::ChooseNucleons not efficient" <<
G4endl;
310 while ( (i <
myA) && (--interationsLeft>0))
317 G4RandFlat::shootArray(jr,prand);
322 aPos.
set((2*arand[jx]-1.), (2*arand[jy]-1.), (2*arand[--jr]-1.));
323 }
while (aPos.
mag2() > 1. );
329 std::vector<G4ThreeVector>::iterator iplace;
330 for( iplace=
places.begin(); iplace!=
places.end() && freeplace;++iplace)
332 delta = *iplace - aPos;
333 freeplace= delta.
mag2() > nd2;
343 G4double eFermi= std::sqrt(
sqr(pFermi) +
sqr(nucMass) ) - nucMass;
354 if (interationsLeft<=0) {
356 "Problem to place nucleons");
373 G4int loopCounterLeft = 10000;
374 for(
G4int ii=1; ii<4; ii++)
382 for(
G4int jj=0; jj < ii; jj++)
384 if( (
theNucleons[ii].GetPosition() -
theNucleons[jj].GetPosition()).mag2() <= nd2 ) {Continue =
true;
break;}
386 }
while( Continue && --loopCounterLeft > 0 );
388 if ( loopCounterLeft <= 0 ) {
390 "Unable to find a good position for the first alpha cluster");
392 loopCounterLeft = 10000;
393 for(
G4int ii=4; ii<8; ii++)
401 for(
G4int jj=0; jj < ii; jj++)
403 if( (
theNucleons[ii].GetPosition() -
theNucleons[jj].GetPosition()).mag2() <= nd2 ) {Continue =
true;
break;}
405 }
while( Continue && --loopCounterLeft > 0 );
407 if ( loopCounterLeft <= 0 ) {
409 "Unable to find a good position for the second alpha cluster");
411 loopCounterLeft = 10000;
412 for(
G4int ii=8; ii<12; ii++)
420 for(
G4int jj=0; jj < ii; jj++)
422 if( (
theNucleons[ii].GetPosition() -
theNucleons[jj].GetPosition()).mag2() <= nd2 ) {Continue =
true;
break;}
424 }
while( Continue && --loopCounterLeft > 0 );
426 if ( loopCounterLeft <= 0 ) {
428 "Unable to find a good position for the third alpha cluster");
454 for (
G4int ntry=0; ntry<1 ; ntry ++ )
456 for (i=0; i <
myA; i++ )
465 if ( eMax >
theNucleons[i].GetDefinition()->GetPDGMass() )
468 fermiM[i] = std::sqrt(pmax2);
469 while ( mom.
mag2() > pmax2 )
478 ed <<
"Nucleus Z A " <<
myZ <<
" " << myA <<
G4endl;
479 ed <<
"proton with eMax=" << eMax <<
G4endl;
480 G4Exception(
"G4Fancy3DNucleus::ChooseFermiMomenta(): difficulty finding proton momentum, set it to (0,0,0)",
499 for ( i=0; i<
myA ; i++ )
501 energy =
theNucleons[i].GetParticleType()->GetPDGMass()
521 if ( sum.
mag() <= PFermi )
533 for (
G4int aNucleon=0; aNucleon < myA-1; aNucleon++) {
534 delta = 2.*((
momentum[aNucleon]*testDir)*testDir);
536 testSums[aNucleon].Fill(delta, delta.
mag(), aNucleon);
543 while ( (sum-
testSums[--index].Vector).mag()>PFermi && index>0)
546 if ( sum.
mag() > (sum-
testSums[index].Vector).mag() ) {
548 sum-=testSums[index].Vector;
552 if ( (sum-
testSums[index].Vector).mag() <= PFermi )
556 for (
G4int aNucleon=0; aNucleon<=index; aNucleon++)
561 && std::abs(
momentum[myA-1].mag() - pTry ) < pBest )
563 pBest=std::abs(
momentum[myA-1].mag() - pTry );
569 G4String text =
"G4Fancy3DNucleus.cc: Logic error in ReduceSum()";
583 if (
fermiM[++swapit] > PFermi )
break;
585 if (swapit == myA-1 )
return false;
void set(double x, double y, double z)
void Init(G4int anA, G4int aZ)
G4VNuclearDensity * theDensity
CLHEP::Hep3Vector G4ThreeVector
std::ostringstream G4ExceptionDescription
G4double Z13(G4int Z) const
static constexpr double MeV
G4double GetFermiMomentum(G4double density)
G4Nucleon * GetNextNucleon()
G4double GetDensity(const G4ThreeVector &aPosition) const
static G4Proton * Proton()
void ChooseFermiMomenta()
G4double GetPDGMass() const
virtual G4double GetRelativeDensity(const G4ThreeVector &aPosition) const =0
const G4ThreeVector & GetPosition() const
static G4double GetBindingEnergy(const G4int A, const G4int Z)
static constexpr double fermi
static G4Pow * GetInstance()
std::vector< G4double > fermiM
HepLorentzRotation & rotateZ(double delta)
virtual G4double GetRadius(const G4double maxRelativeDenisty) const =0
void DoLorentzContraction(const G4LorentzVector &theBoost)
std::vector< G4Fancy3DNucleusHelper > testSums
G4double GetOuterRadius()
G4ThreeVector GetMomentum(G4double density, G4double maxMomentum=-1.)
std::vector< G4Nucleon > theNucleons
void Init(G4int theA, G4int theZ)
const std::vector< G4Nucleon > & GetNucleons()
ThreeVector shoot(const G4int Ap, const G4int Af)
bool G4Fancy3DNucleusHelperForSortInZ(const G4Nucleon &nuc1, const G4Nucleon &nuc2)
std::vector< G4ThreeVector > momentum
G4double CoulombBarrier()
const G4VNuclearDensity * GetNuclearDensity() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double GetNuclearRadius()
G4double excitationEnergy
static G4Neutron * Neutron()
G4GLOB_DLL std::ostream G4cout
static constexpr double pi
CLHEP::HepLorentzVector G4LorentzVector
HepLorentzRotation & rotateY(double delta)
std::vector< G4ThreeVector > places
void DoLorentzBoost(const G4LorentzVector &theBoost)
static constexpr double GeV
void DoTranslation(const G4ThreeVector &theShift)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments