Geant4  10.03
G4Fancy3DNucleus.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // ------------------------------------------------------------
27 // GEANT 4 class implementation file
28 //
29 // ---------------- G4Fancy3DNucleus ----------------
30 // by Gunter Folger, May 1998.
31 // class for a 3D nucleus, arranging nucleons in space and momentum.
32 // ------------------------------------------------------------
33 // 20110805 M. Kelsey -- Remove C-style array (pointer) of G4Nucleons,
34 // make vector a container of objects. Move Helper class
35 // to .hh. Move testSums, places, momentum and fermiM to
36 // class data members for reuse.
37 
38 #include <algorithm>
39 
40 #include "G4Fancy3DNucleus.hh"
42 #include "G4NuclearFermiDensity.hh"
44 #include "G4NucleiProperties.hh"
45 #include "G4Nucleon.hh"
46 #include "G4SystemOfUnits.hh"
47 #include "Randomize.hh"
48 #include "G4ios.hh"
49 #include "G4Pow.hh"
50 #include "G4HadronicException.hh"
51 
52 
54  : myA(0), myZ(0), theNucleons(250), currentNucleon(-1), theDensity(0),
55  nucleondistance(0.8*fermi),excitationEnergy(0.),
56  places(250), momentum(250), fermiM(250), testSums(250)
57 {
58 //G4cout <<"G4Fancy3DNucleus::G4Fancy3DNucleus()"<<G4endl;
59 }
60 
62 {
63  if(theDensity) delete theDensity;
64 }
65 
66 #if defined(NON_INTEGER_A_Z)
68 {
69  G4int intZ = G4int(theZ);
70  G4int intA= ( G4UniformRand()>theA-G4int(theA) ) ? G4int(theA) : G4int(theA)+1;
71  // forward to integer Init()
72  Init(intA, intZ);
73 
74 }
75 #endif
76 
78 {
79 // G4cout << "G4Fancy3DNucleus::Init(theA, theZ) called"<<G4endl;
80  currentNucleon=-1;
81  theNucleons.clear();
82 
83  myZ = theZ;
84  myA= theA;
86 
87  theNucleons.resize(myA); // Pre-loads vector with empty elements
88 
89 // G4cout << "myA, myZ" << myA << ", " << myZ << G4endl;
90 
91  if(theDensity) delete theDensity;
92  if ( myA < 17 ) {
94  } else {
96  }
97 
98  theFermi.Init(myA, myZ);
99 
100  ChooseNucleons();
101 
102  ChoosePositions();
103 
104 // CenterNucleons(); // This would introduce a bias
105 
107 
108  G4double Ebinding= BindingEnergy()/myA;
109 
110  for (G4int aNucleon=0; aNucleon < myA; aNucleon++)
111  {
112  theNucleons[aNucleon].SetBindingEnergy(Ebinding);
113  }
114 
115 
116  return;
117 }
118 
120 {
121  currentNucleon=0;
122  return (theNucleons.size()>0);
123 }
124 
125 // Returns by pointer; null pointer indicates end of loop
127 {
128  return ( (currentNucleon>=0 && currentNucleon<myA) ?
129  &theNucleons[currentNucleon++] : 0 );
130 }
131 
132 const std::vector<G4Nucleon> & G4Fancy3DNucleus::GetNucleons()
133 {
134  return theNucleons;
135 }
136 
137 
138 // Class-scope function to sort nucleons by Z coordinate
140 {
141  return nuc1.GetPosition().z() < nuc2.GetPosition().z();
142 }
143 
144 void G4Fancy3DNucleus::SortNucleonsIncZ() // on increased Z-coordinates Uzhi 29.08.08
145 {
146  if (theNucleons.size() < 2 ) return; // Avoid unnecesary work
147 
148  std::sort(theNucleons.begin(), theNucleons.end(),
150 }
151 
152 void G4Fancy3DNucleus::SortNucleonsDecZ() // on decreased Z-coordinates Uzhi 29.08.08
153 {
154  if (theNucleons.size() < 2 ) return; // Avoid unnecessary work
156 
157  std::reverse(theNucleons.begin(), theNucleons.end());
158 }
159 
160 
162 {
164 }
165 
166 
168 {
169  return GetNuclearRadius(0.5);
170 }
171 
173 {
174  return theDensity->GetRadius(maxRelativeDensity);
175 }
176 
178 {
179  G4double maxradius2=0;
180 
181  for (int i=0; i<myA; i++)
182  {
183  if ( theNucleons[i].GetPosition().mag2() > maxradius2 )
184  {
185  maxradius2=theNucleons[i].GetPosition().mag2();
186  }
187  }
188  return std::sqrt(maxradius2)+nucleondistance;
189 }
190 
192 {
193  return myZ*G4Proton::Proton()->GetPDGMass() +
195  BindingEnergy();
196 }
197 
198 
199 
201 {
202  for (G4int i=0; i<myA; i++){
203  theNucleons[i].Boost(theBoost);
204  }
205 }
206 
208 {
209  for (G4int i=0; i<myA; i++){
210  theNucleons[i].Boost(theBeta);
211  }
212 }
213 
215 {
216  G4double beta2=theBeta.mag2();
217  if (beta2 > 0) {
218  G4double factor=(1-std::sqrt(1-beta2))/beta2; // (gamma-1)/gamma/beta**2
219  G4ThreeVector rprime;
220  for (G4int i=0; i< myA; i++) {
221  rprime = theNucleons[i].GetPosition() -
222  factor * (theBeta*theNucleons[i].GetPosition()) * theBeta;
223  theNucleons[i].SetPosition(rprime);
224  }
225  }
226 }
227 
229 {
230  if (theBoost.e() !=0 ) {
231  G4ThreeVector beta = theBoost.vect()/theBoost.e();
232  DoLorentzContraction(beta);
233  }
234 }
235 
236 
237 
239 {
240  G4ThreeVector center;
241 
242  for (G4int i=0; i<myA; i++ )
243  {
244  center+=theNucleons[i].GetPosition();
245  }
246  center /= -myA;
247  DoTranslation(center);
248 }
249 
251 {
252  G4ThreeVector tempV;
253  for (G4int i=0; i<myA; i++ )
254  {
255  tempV = theNucleons[i].GetPosition() + theShift;
256  theNucleons[i].SetPosition(tempV);
257  }
258 }
259 
261 {
262  return theDensity;
263 }
264 
265 //----------------------- private Implementation Methods-------------
266 
268 {
269  G4int protons=0,nucleons=0;
270 
271  while (nucleons < myA ) /* Loop checking, 30-Oct-2015, G.Folger */
272  {
273  if ( protons < myZ && G4UniformRand() < (G4double)(myZ-protons)/(G4double)(myA-nucleons) )
274  {
275  protons++;
276  theNucleons[nucleons++].SetParticleType(G4Proton::Proton());
277  }
278  else if ( (nucleons-protons) < (myA-myZ) )
279  {
280  theNucleons[nucleons++].SetParticleType(G4Neutron::Neutron());
281  }
282  else G4cout << "G4Fancy3DNucleus::ChooseNucleons not efficient" << G4endl;
283  }
284  return;
285 }
286 
288 {
289  G4int i=0;
290  G4ThreeVector aPos, delta;
291  G4bool freeplace;
292  const G4double nd2=sqr(nucleondistance);
293  G4double maxR=GetNuclearRadius(0.001); // there are no nucleons at a
294  // relative Density of 0.01
295  G4int jr=0;
296  G4int jx,jy;
297  G4double arand[600];
298  G4double *prand=arand;
299 
300  places.clear(); // Reset data buffer
301  G4int interationsLeft=1000*myA;
302  while ( (i < myA) && (--interationsLeft>0)) /* Loop checking, 30-Oct-2015, G.Folger */
303  {
304  do
305  {
306  if ( jr < 3 )
307  {
308  jr=std::min(600,9*(myA - i));
309  G4RandFlat::shootArray(jr,prand);
310  //CLHEP::RandFlat::shootArray(jr, prand );
311  }
312  jx=--jr;
313  jy=--jr;
314  aPos.set((2*arand[jx]-1.), (2*arand[jy]-1.), (2*arand[--jr]-1.));
315  } while (aPos.mag2() > 1. ); /* Loop checking, 30-Oct-2015, G.Folger */
316  aPos *=maxR;
317  G4double density=theDensity->GetRelativeDensity(aPos);
318  if (G4UniformRand() < density)
319  {
320  freeplace= true;
321  std::vector<G4ThreeVector>::iterator iplace;
322  for( iplace=places.begin(); iplace!=places.end() && freeplace;++iplace)
323  {
324  delta = *iplace - aPos;
325  freeplace= delta.mag2() > nd2;
326  }
327 
328  if ( freeplace )
329  {
331  // protons must at least have binding energy of CoulombBarrier, so
332  // assuming the Fermi energy corresponds to a potential, we must place these such
333  // that the Fermi Energy > CoulombBarrier
334  if (theNucleons[i].GetDefinition() == G4Proton::Proton())
335  {
336  G4double nucMass = theNucleons[i].GetDefinition()->GetPDGMass();
337  G4double eFermi= std::sqrt( sqr(pFermi) + sqr(nucMass) )
338  - nucMass;
339  if (eFermi <= CoulombBarrier() ) freeplace=false;
340  }
341  }
342  if ( freeplace )
343  {
344  theNucleons[i].SetPosition(aPos);
345  places.push_back(aPos);
346  ++i;
347  }
348  }
349  }
350  if (interationsLeft<=0) {
351  G4Exception("model/util/G4Fancy3DNucleus.cc", "mod_util001", FatalException,
352  "Problem to place nucleons");
353  }
354 
355 }
356 
358 {
359  G4int i;
360  G4double density;
361 
362  // Pre-allocate buffers for filling by index
363  momentum.resize(myA, G4ThreeVector(0.,0.,0.));
364  fermiM.resize(myA, 0.*GeV);
365 
366  for (G4int ntry=0; ntry<1 ; ntry ++ )
367  {
368  for (i=0; i < myA; i++ ) // momenta for all, including last, in case we swap nucleons
369  {
370  density = theDensity->GetDensity(theNucleons[i].GetPosition());
371  fermiM[i] = theFermi.GetFermiMomentum(density);
372  G4ThreeVector mom=theFermi.GetMomentum(density);
373  if (theNucleons[i].GetDefinition() == G4Proton::Proton())
374  {
375  G4double eMax = std::sqrt(sqr(fermiM[i]) +sqr(theNucleons[i].GetDefinition()->GetPDGMass()) )
376  - CoulombBarrier();
377  if ( eMax > theNucleons[i].GetDefinition()->GetPDGMass() )
378  {
379  G4double pmax2= sqr(eMax) - sqr(theNucleons[i].GetDefinition()->GetPDGMass());
380  fermiM[i] = std::sqrt(pmax2);
381  while ( mom.mag2() > pmax2 ) /* Loop checking, 30-Oct-2015, G.Folger */
382  {
383  mom=theFermi.GetMomentum(density, fermiM[i]);
384  }
385  } else
386  {
387  G4cerr << "G4Fancy3DNucleus: difficulty finding proton momentum" << G4endl;
388  mom=G4ThreeVector(0,0,0);
389  }
390 
391  }
392  momentum[i]= mom;
393  }
394 
395  if ( ReduceSum() ) break;
396 // G4cout <<" G4FancyNucleus: iterating to find momenta: "<< ntry<< G4endl;
397  }
398 
399 // G4ThreeVector sum;
400 // for (G4int index=0; index<myA;sum+=momentum[index++])
401 // ;
402 // G4cout << "final sum / mag() " << sum << " / " << sum.mag() << G4endl;
403 
405  for ( i=0; i< myA ; i++ )
406  {
407  energy = theNucleons[i].GetParticleType()->GetPDGMass()
408  - BindingEnergy()/myA;
409  G4LorentzVector tempV(momentum[i],energy);
410  theNucleons[i].SetMomentum(tempV);
411  // GF 11-05-2011: set BindingEnergy to be T of Nucleon with p , ~ p**2/2m
412  //theNucleons[i].SetBindingEnergy(
413  // 0.5*sqr(fermiM[i])/theNucleons[i].GetParticleType()->GetPDGMass());
414  }
415 }
416 
417 
419 {
420  G4ThreeVector sum;
421  G4double PFermi=fermiM[myA-1];
422 
423  for (G4int i=0; i < myA-1 ; i++ )
424  { sum+=momentum[i]; }
425 
426 // check if have to do anything at all..
427  if ( sum.mag() <= PFermi )
428  {
429  momentum[myA-1]=-sum;
430  return true;
431  }
432 
433 // find all possible changes in momentum, changing only the component parallel to sum
434  G4ThreeVector testDir=sum.unit();
435  testSums.clear();
436  testSums.resize(myA-1); // Allocate block for filling below
437 
438  G4ThreeVector delta;
439  for (G4int aNucleon=0; aNucleon < myA-1; aNucleon++) {
440  delta = 2.*((momentum[aNucleon]*testDir)*testDir);
441 
442  testSums[aNucleon].Fill(delta, delta.mag(), aNucleon);
443  }
444 
445  std::sort(testSums.begin(), testSums.end());
446 
447 // reduce Momentum Sum until the next would be allowed.
448  G4int index=testSums.size();
449  while ( (sum-testSums[--index].Vector).mag()>PFermi && index>0) /* Loop checking, 30-Oct-2015, G.Folger */
450  {
451  // Only take one which improve, ie. don't change sign and overshoot...
452  if ( sum.mag() > (sum-testSums[index].Vector).mag() ) {
453  momentum[testSums[index].Index]-=testSums[index].Vector;
454  sum-=testSums[index].Vector;
455  }
456  }
457 
458  if ( (sum-testSums[index].Vector).mag() <= PFermi )
459  {
460  G4int best=-1;
461  G4double pBest=2*PFermi; // anything larger than PFermi
462  for ( G4int aNucleon=0; aNucleon<=index; aNucleon++)
463  {
464  // find the momentum closest to choosen momentum for last Nucleon.
465  G4double pTry=(testSums[aNucleon].Vector-sum).mag();
466  if ( pTry < PFermi
467  && std::abs(momentum[myA-1].mag() - pTry ) < pBest )
468  {
469  pBest=std::abs(momentum[myA-1].mag() - pTry );
470  best=aNucleon;
471  }
472  }
473  if ( best < 0 )
474  {
475  G4String text = "G4Fancy3DNucleus.cc: Logic error in ReduceSum()";
476  throw G4HadronicException(__FILE__, __LINE__, text);
477  }
478  momentum[testSums[best].Index]-=testSums[best].Vector;
479  momentum[myA-1]=testSums[best].Vector-sum;
480 
481  return true;
482 
483  }
484 
485  // try to compensate momentum using another Nucleon....
486  G4int swapit=-1;
487  while (swapit<myA-1) /* Loop checking, 30-Oct-2015, G.Folger */
488  {
489  if ( fermiM[++swapit] > PFermi ) break;
490  }
491  if (swapit == myA-1 ) return false;
492 
493  // Now we have a nucleon with a bigger Fermi Momentum.
494  // Exchange with last nucleon.. and iterate.
495  std::swap(theNucleons[swapit], theNucleons[myA-1]);
496  std::swap(momentum[swapit], momentum[myA-1]);
497  std::swap(fermiM[swapit], fermiM[myA-1]);
498  return ReduceSum();
499 }
500 
502 {
503  static const G4double cfactor = (1.44/1.14) * MeV;
504  return cfactor*myZ/(1.0 + G4Pow::GetInstance()->Z13(myA));
505 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
std::vector< G4ThreeVector > momentum
CLHEP::Hep3Vector G4ThreeVector
G4double BindingEnergy()
G4double CoulombBarrier()
void DoTranslation(const G4ThreeVector &theShift)
std::vector< G4double > fermiM
void DoLorentzBoost(const G4LorentzVector &theBoost)
virtual const G4ThreeVector & GetPosition() const
Definition: G4Nucleon.hh:68
std::vector< G4Fancy3DNucleusHelper > testSums
int G4int
Definition: G4Types.hh:78
std::vector< G4ThreeVector > places
virtual G4double GetRadius(const G4double maxRelativeDenisty) const =0
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4double Z13(G4int Z) const
Definition: G4Pow.hh:127
std::vector< G4Nucleon > theNucleons
bool G4bool
Definition: G4Types.hh:79
bool G4Fancy3DNucleusHelperForSortInZ(const G4Nucleon &nuc1, const G4Nucleon &nuc2)
void Init(G4int theA, G4int theZ)
virtual G4double GetRelativeDensity(const G4ThreeVector &aPosition) const =0
G4double GetFermiMomentum(G4double density)
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4double GetNuclearRadius()
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
const std::vector< G4Nucleon > & GetNucleons()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void Init(G4int anA, G4int aZ)
G4double GetPDGMass() const
G4double GetDensity(const G4ThreeVector &aPosition) const
static G4double GetBindingEnergy(const G4int A, const G4int Z)
G4double energy(const ThreeVector &p, const G4double m)
void DoLorentzContraction(const G4LorentzVector &theBoost)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static constexpr double GeV
Definition: G4SIunits.hh:217
G4FermiMomentum theFermi
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
G4ThreeVector GetMomentum(G4double density, G4double maxMomentum=-1.)
G4double GetOuterRadius()
T sqr(const T &x)
Definition: templates.hh:145
const G4VNuclearDensity * GetNuclearDensity() const
double G4double
Definition: G4Types.hh:76
const G4double nucleondistance
static constexpr double fermi
Definition: G4SIunits.hh:103
G4Nucleon * GetNextNucleon()
G4VNuclearDensity * theDensity
G4GLOB_DLL std::ostream G4cerr
CLHEP::HepLorentzVector G4LorentzVector