Geant4_10
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 "G4HadronicException.hh"
50 
51 
53  : myA(0), myZ(0), theNucleons(250), currentNucleon(-1), theDensity(0),
54  nucleondistance(0.8*fermi),excitationEnergy(0.),
55  places(250), momentum(250), fermiM(250), testSums(250)
56 {
57 //G4cout <<"G4Fancy3DNucleus::G4Fancy3DNucleus()"<<G4endl;
58 }
59 
61 {
62  if(theDensity) delete theDensity;
63 }
64 
65 #if defined(NON_INTEGER_A_Z)
67 {
68  G4int intZ = G4int(theZ);
69  G4int intA= ( G4UniformRand()>theA-G4int(theA) ) ? G4int(theA) : G4int(theA)+1;
70  // forward to integer Init()
71  Init(intA, intZ);
72 
73 }
74 #endif
75 
77 {
78 // G4cout << "G4Fancy3DNucleus::Init(theA, theZ) called"<<G4endl;
79  currentNucleon=-1;
80  theNucleons.clear();
81 
82  myZ = theZ;
83  myA= theA;
84  excitationEnergy=0;
85 
86  theNucleons.resize(myA); // Pre-loads vector with empty elements
87 
88 // G4cout << "myA, myZ" << myA << ", " << myZ << G4endl;
89 
90  if(theDensity) delete theDensity;
91  if ( myA < 17 ) {
92  theDensity = new G4NuclearShellModelDensity(myA, myZ);
93  } else {
94  theDensity = new G4NuclearFermiDensity(myA, myZ);
95  }
96 
97  theFermi.Init(myA, myZ);
98 
99  ChooseNucleons();
100 
101  ChoosePositions();
102 
103 // CenterNucleons(); // This would introduce a bias
104 
105  ChooseFermiMomenta();
106 
107  G4double Ebinding= BindingEnergy()/myA;
108 
109  for (G4int aNucleon=0; aNucleon < myA; aNucleon++)
110  {
111  theNucleons[aNucleon].SetBindingEnergy(Ebinding);
112  }
113 
114 
115  return;
116 }
117 
119 {
120  currentNucleon=0;
121  return (theNucleons.size()>0);
122 }
123 
124 // Returns by pointer; null pointer indicates end of loop
126 {
127  return ( (currentNucleon>=0 && currentNucleon<myA) ?
128  &theNucleons[currentNucleon++] : 0 );
129 }
130 
131 const std::vector<G4Nucleon> & G4Fancy3DNucleus::GetNucleons()
132 {
133  return theNucleons;
134 }
135 
136 
137 // Class-scope function to sort nucleons by Z coordinate
139 {
140  return nuc1.GetPosition().z() < nuc2.GetPosition().z();
141 }
142 
143 void G4Fancy3DNucleus::SortNucleonsIncZ() // on increased Z-coordinates Uzhi 29.08.08
144 {
145  if (theNucleons.size() < 2 ) return; // Avoid unnecesary work
146 
147  std::sort(theNucleons.begin(), theNucleons.end(),
149 }
150 
151 void G4Fancy3DNucleus::SortNucleonsDecZ() // on decreased Z-coordinates Uzhi 29.08.08
152 {
153  if (theNucleons.size() < 2 ) return; // Avoid unnecessary work
155 
156  std::reverse(theNucleons.begin(), theNucleons.end());
157 }
158 
159 
160 G4double G4Fancy3DNucleus::BindingEnergy()
161 {
162  return G4NucleiProperties::GetBindingEnergy(myA,myZ);
163 }
164 
165 
167 {
168  return GetNuclearRadius(0.5);
169 }
170 
172 {
173  return theDensity->GetRadius(maxRelativeDensity);
174 }
175 
177 {
178  G4double maxradius2=0;
179 
180  for (int i=0; i<myA; i++)
181  {
182  if ( theNucleons[i].GetPosition().mag2() > maxradius2 )
183  {
184  maxradius2=theNucleons[i].GetPosition().mag2();
185  }
186  }
187  return std::sqrt(maxradius2)+nucleondistance;
188 }
189 
191 {
192  return myZ*G4Proton::Proton()->GetPDGMass() +
193  (myA-myZ)*G4Neutron::Neutron()->GetPDGMass() -
194  BindingEnergy();
195 }
196 
197 
198 
200 {
201  for (G4int i=0; i<myA; i++){
202  theNucleons[i].Boost(theBoost);
203  }
204 }
205 
207 {
208  for (G4int i=0; i<myA; i++){
209  theNucleons[i].Boost(theBeta);
210  }
211 }
212 
214 {
215  G4double beta2=theBeta.mag2();
216  if (beta2 > 0) {
217  G4double factor=(1-std::sqrt(1-beta2))/beta2; // (gamma-1)/gamma/beta**2
218  G4ThreeVector rprime;
219  for (G4int i=0; i< myA; i++) {
220  rprime = theNucleons[i].GetPosition() -
221  factor * (theBeta*theNucleons[i].GetPosition()) * theBeta;
222  theNucleons[i].SetPosition(rprime);
223  }
224  }
225 }
226 
228 {
229  if (theBoost.e() !=0 ) {
230  G4ThreeVector beta = theBoost.vect()/theBoost.e();
231  DoLorentzContraction(beta);
232  }
233 }
234 
235 
236 
238 {
239  G4ThreeVector center;
240 
241  for (G4int i=0; i<myA; i++ )
242  {
243  center+=theNucleons[i].GetPosition();
244  }
245  center /= -myA;
246  DoTranslation(center);
247 }
248 
250 {
251  G4ThreeVector tempV;
252  for (G4int i=0; i<myA; i++ )
253  {
254  tempV = theNucleons[i].GetPosition() + theShift;
255  theNucleons[i].SetPosition(tempV);
256  }
257 }
258 
260 {
261  return theDensity;
262 }
263 
264 //----------------------- private Implementation Methods-------------
265 
266 void G4Fancy3DNucleus::ChooseNucleons()
267 {
268  G4int protons=0,nucleons=0;
269 
270  while (nucleons < myA )
271  {
272  if ( protons < myZ && G4UniformRand() < (G4double)(myZ-protons)/(G4double)(myA-nucleons) )
273  {
274  protons++;
275  theNucleons[nucleons++].SetParticleType(G4Proton::Proton());
276  }
277  else if ( (nucleons-protons) < (myA-myZ) )
278  {
279  theNucleons[nucleons++].SetParticleType(G4Neutron::Neutron());
280  }
281  else G4cout << "G4Fancy3DNucleus::ChooseNucleons not efficient" << G4endl;
282  }
283  return;
284 }
285 
286 void G4Fancy3DNucleus::ChoosePositions()
287 {
288  G4int i=0;
289  G4ThreeVector aPos, delta;
290  G4bool freeplace;
291  static G4ThreadLocal G4double *nd2_G4MT_TLS_ = 0 ; if (!nd2_G4MT_TLS_) {nd2_G4MT_TLS_ = new G4double ; *nd2_G4MT_TLS_= sqr(nucleondistance) ; } G4double &nd2 = *nd2_G4MT_TLS_;
292  G4double maxR=GetNuclearRadius(0.001); // there are no nucleons at a
293  // relative Density of 0.01
294  G4int jr=0;
295  G4int jx,jy;
296  G4double arand[600];
297  G4double *prand=arand;
298 
299  places.clear(); // Reset data buffer
300 
301  while ( i < myA )
302  {
303  do
304  {
305  if ( jr < 3 )
306  {
307  jr=std::min(600,9*(myA - i));
308  G4RandFlat::shootArray(jr,prand);
309  //CLHEP::RandFlat::shootArray(jr, prand );
310  }
311  jx=--jr;
312  jy=--jr;
313  aPos.set((2*arand[jx]-1.), (2*arand[jy]-1.), (2*arand[--jr]-1.));
314  } while (aPos.mag2() > 1. );
315  aPos *=maxR;
316  G4double density=theDensity->GetRelativeDensity(aPos);
317  if (G4UniformRand() < density)
318  {
319  freeplace= true;
320  std::vector<G4ThreeVector>::iterator iplace;
321  for( iplace=places.begin(); iplace!=places.end() && freeplace;++iplace)
322  {
323  delta = *iplace - aPos;
324  freeplace= delta.mag2() > nd2;
325  }
326 
327  if ( freeplace )
328  {
329  G4double pFermi=theFermi.GetFermiMomentum(theDensity->GetDensity(aPos));
330  // protons must at least have binding energy of CoulombBarrier, so
331  // assuming the Fermi energy corresponds to a potential, we must place these such
332  // that the Fermi Energy > CoulombBarrier
333  if (theNucleons[i].GetDefinition() == G4Proton::Proton())
334  {
335  G4double nucMass = theNucleons[i].GetDefinition()->GetPDGMass();
336  G4double eFermi= std::sqrt( sqr(pFermi) + sqr(nucMass) )
337  - nucMass;
338  if (eFermi <= CoulombBarrier() ) freeplace=false;
339  }
340  }
341  if ( freeplace )
342  {
343  theNucleons[i].SetPosition(aPos);
344  places.push_back(aPos);
345  ++i;
346  }
347  }
348  }
349 
350 }
351 
352 void G4Fancy3DNucleus::ChooseFermiMomenta()
353 {
354  G4int i;
356 
357  // Pre-allocate buffers for filling by index
358  momentum.resize(myA, G4ThreeVector(0.,0.,0.));
359  fermiM.resize(myA, 0.*GeV);
360 
361  for (G4int ntry=0; ntry<1 ; ntry ++ )
362  {
363  for (i=0; i < myA; i++ ) // momenta for all, including last, in case we swap nucleons
364  {
365  density = theDensity->GetDensity(theNucleons[i].GetPosition());
366  fermiM[i] = theFermi.GetFermiMomentum(density);
367  G4ThreeVector mom=theFermi.GetMomentum(density);
368  if (theNucleons[i].GetDefinition() == G4Proton::Proton())
369  {
370  G4double eMax = std::sqrt(sqr(fermiM[i]) +sqr(theNucleons[i].GetDefinition()->GetPDGMass()) )
371  - CoulombBarrier();
372  if ( eMax > theNucleons[i].GetDefinition()->GetPDGMass() )
373  {
374  G4double pmax2= sqr(eMax) - sqr(theNucleons[i].GetDefinition()->GetPDGMass());
375  fermiM[i] = std::sqrt(pmax2);
376  while ( mom.mag2() > pmax2 )
377  {
378  mom=theFermi.GetMomentum(density, fermiM[i]);
379  }
380  } else
381  {
382  G4cerr << "G4Fancy3DNucleus: difficulty finding proton momentum" << G4endl;
383  mom=G4ThreeVector(0,0,0);
384  }
385 
386  }
387  momentum[i]= mom;
388  }
389 
390  if ( ReduceSum() ) break;
391 // G4cout <<" G4FancyNucleus: iterating to find momenta: "<< ntry<< G4endl;
392  }
393 
394 // G4ThreeVector sum;
395 // for (G4int index=0; index<myA;sum+=momentum[index++])
396 // ;
397 // G4cout << "final sum / mag() " << sum << " / " << sum.mag() << G4endl;
398 
400  for ( i=0; i< myA ; i++ )
401  {
402  energy = theNucleons[i].GetParticleType()->GetPDGMass()
403  - BindingEnergy()/myA;
404  G4LorentzVector tempV(momentum[i],energy);
405  theNucleons[i].SetMomentum(tempV);
406  // GF 11-05-2011: set BindingEnergy to be T of Nucleon with p , ~ p**2/2m
407  //theNucleons[i].SetBindingEnergy(
408  // 0.5*sqr(fermiM[i])/theNucleons[i].GetParticleType()->GetPDGMass());
409  }
410 }
411 
412 
413 G4bool G4Fancy3DNucleus::ReduceSum()
414 {
415  G4ThreeVector sum;
416  G4double PFermi=fermiM[myA-1];
417 
418  for (G4int i=0; i < myA-1 ; i++ )
419  { sum+=momentum[i]; }
420 
421 // check if have to do anything at all..
422  if ( sum.mag() <= PFermi )
423  {
424  momentum[myA-1]=-sum;
425  return true;
426  }
427 
428 // find all possible changes in momentum, changing only the component parallel to sum
429  G4ThreeVector testDir=sum.unit();
430  testSums.clear();
431  testSums.resize(myA-1); // Allocate block for filling below
432 
433  G4ThreeVector delta;
434  for (G4int aNucleon=0; aNucleon < myA-1; aNucleon++) {
435  delta = 2.*((momentum[aNucleon]*testDir)*testDir);
436 
437  testSums[aNucleon].Fill(delta, delta.mag(), aNucleon);
438  }
439 
440  std::sort(testSums.begin(), testSums.end());
441 
442 // reduce Momentum Sum until the next would be allowed.
443  G4int index=testSums.size();
444  while ( (sum-testSums[--index].Vector).mag()>PFermi && index>0)
445  {
446  // Only take one which improve, ie. don't change sign and overshoot...
447  if ( sum.mag() > (sum-testSums[index].Vector).mag() ) {
448  momentum[testSums[index].Index]-=testSums[index].Vector;
449  sum-=testSums[index].Vector;
450  }
451  }
452 
453  if ( (sum-testSums[index].Vector).mag() <= PFermi )
454  {
455  G4int best=-1;
456  G4double pBest=2*PFermi; // anything larger than PFermi
457  for ( G4int aNucleon=0; aNucleon<=index; aNucleon++)
458  {
459  // find the momentum closest to choosen momentum for last Nucleon.
460  G4double pTry=(testSums[aNucleon].Vector-sum).mag();
461  if ( pTry < PFermi
462  && std::abs(momentum[myA-1].mag() - pTry ) < pBest )
463  {
464  pBest=std::abs(momentum[myA-1].mag() - pTry );
465  best=aNucleon;
466  }
467  }
468  if ( best < 0 )
469  {
470  G4String text = "G4Fancy3DNucleus.cc: Logic error in ReduceSum()";
471  throw G4HadronicException(__FILE__, __LINE__, text);
472  }
473  momentum[testSums[best].Index]-=testSums[best].Vector;
474  momentum[myA-1]=testSums[best].Vector-sum;
475 
476  return true;
477 
478  }
479 
480  // try to compensate momentum using another Nucleon....
481  G4int swapit=-1;
482  while (swapit<myA-1)
483  {
484  if ( fermiM[++swapit] > PFermi ) break;
485  }
486  if (swapit == myA-1 ) return false;
487 
488  // Now we have a nucleon with a bigger Fermi Momentum.
489  // Exchange with last nucleon.. and iterate.
490  std::swap(theNucleons[swapit], theNucleons[myA-1]);
491  std::swap(momentum[swapit], momentum[myA-1]);
492  std::swap(fermiM[swapit], fermiM[myA-1]);
493  return ReduceSum();
494 }
495 
497 {
498  G4double coulombBarrier = (1.44/1.14) * MeV * myZ / (1.0 + std::pow(G4double(myA),1./3.));
499  return coulombBarrier;
500 }
void set(double x, double y, double z)
CLHEP::Hep3Vector G4ThreeVector
Int_t index
Definition: macro.C:9
G4double CoulombBarrier()
void DoTranslation(const G4ThreeVector &theShift)
subroutine sort(A, N)
Definition: dpm25nuc7.f:4670
void DoLorentzBoost(const G4LorentzVector &theBoost)
virtual const G4ThreeVector & GetPosition() const
Definition: G4Nucleon.hh:68
#define G4ThreadLocal
Definition: tls.hh:52
int G4int
Definition: G4Types.hh:78
double z() const
virtual G4double GetRadius(const G4double maxRelativeDenisty) const =0
G4double density
Definition: TRTMaterials.hh:39
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
double energy
Definition: plottest35.C:25
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)
void swap(shared_ptr< P > &, shared_ptr< P > &)
Definition: memory.h:1247
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4double GetNuclearRadius()
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
const std::vector< G4Nucleon > & GetNucleons()
void Init(G4int anA, G4int aZ)
G4double GetPDGMass() const
G4double GetDensity(const G4ThreeVector &aPosition) const
static G4double GetBindingEnergy(const G4int A, const G4int Z)
Hep3Vector unit() const
void DoLorentzContraction(const G4LorentzVector &theBoost)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
double mag2() const
#define G4endl
Definition: G4ios.hh:61
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
double mag() const
G4Nucleon * GetNextNucleon()
G4GLOB_DLL std::ostream G4cerr