ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Fancy3DNucleus.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Fancy3DNucleus.cc
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 #include "Randomize.hh"
53 #include "G4ThreeVector.hh"
54 #include "G4RandomDirection.hh"
55 #include "G4LorentzRotation.hh"
56 #include "G4RotationMatrix.hh"
57 #include "G4PhysicalConstants.hh"
58 
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)
63 {
64 }
65 
67 {
68  if(theDensity) delete theDensity;
69 }
70 
71 #if defined(NON_INTEGER_A_Z)
73 {
74  G4int intZ = G4int(theZ);
75  G4int intA= ( G4UniformRand()>theA-G4int(theA) ) ? G4int(theA) : G4int(theA)+1;
76  // forward to integer Init()
77  Init(intA, intZ);
78 
79 }
80 #endif
81 
83 {
84  currentNucleon=-1;
85  theNucleons.clear();
86  nucleondistance = 0.8*fermi;
87  places.clear();
88  momentum.clear();
89  fermiM.clear();
90  testSums.clear();
91 
92  myZ = theZ;
93  myA= theA;
95 
96  theNucleons.resize(myA); // Pre-loads vector with empty elements
97 
98  if(theDensity) delete theDensity;
99  if ( myA < 17 ) {
101  if( myA == 12 ) nucleondistance=0.9*fermi;
102  } else {
104  }
105 
106  theFermi.Init(myA, myZ);
107 
108  ChooseNucleons();
109 
110  ChoosePositions();
111 
112  if( myA == 12 ) CenterNucleons(); // This would introduce a bias
113 
115 
116  G4double Ebinding= BindingEnergy()/myA;
117 
118  for (G4int aNucleon=0; aNucleon < myA; aNucleon++)
119  {
120  theNucleons[aNucleon].SetBindingEnergy(Ebinding);
121  }
122 
123  return;
124 }
125 
127 {
128  currentNucleon=0;
129  return (theNucleons.size()>0);
130 }
131 
132 // Returns by pointer; null pointer indicates end of loop
134 {
135  return ( (currentNucleon>=0 && currentNucleon<myA) ?
136  &theNucleons[currentNucleon++] : 0 );
137 }
138 
139 const std::vector<G4Nucleon> & G4Fancy3DNucleus::GetNucleons()
140 {
141  return theNucleons;
142 }
143 
144 
145 // Class-scope function to sort nucleons by Z coordinate
147 {
148  return nuc1.GetPosition().z() < nuc2.GetPosition().z();
149 }
150 
152 {
153  if (theNucleons.size() < 2 ) return; // Avoid unnecesary work
154 
155  std::sort(theNucleons.begin(), theNucleons.end(),
157 }
158 
160 {
161  if (theNucleons.size() < 2 ) return; // Avoid unnecessary work
163 
164  std::reverse(theNucleons.begin(), theNucleons.end());
165 }
166 
167 
169 {
171 }
172 
173 
175 {
176  return GetNuclearRadius(0.5);
177 }
178 
180 {
181  return theDensity->GetRadius(maxRelativeDensity);
182 }
183 
185 {
186  G4double maxradius2=0;
187 
188  for (int i=0; i<myA; i++)
189  {
190  if ( theNucleons[i].GetPosition().mag2() > maxradius2 )
191  {
192  maxradius2=theNucleons[i].GetPosition().mag2();
193  }
194  }
195  return std::sqrt(maxradius2)+nucleondistance;
196 }
197 
199 {
200  return myZ*G4Proton::Proton()->GetPDGMass() +
202  BindingEnergy();
203 }
204 
205 
206 
208 {
209  for (G4int i=0; i<myA; i++){
210  theNucleons[i].Boost(theBoost);
211  }
212 }
213 
215 {
216  for (G4int i=0; i<myA; i++){
217  theNucleons[i].Boost(theBeta);
218  }
219 }
220 
222 {
223  G4double beta2=theBeta.mag2();
224  if (beta2 > 0) {
225  G4double factor=(1-std::sqrt(1-beta2))/beta2; // (gamma-1)/gamma/beta**2
226  G4ThreeVector rprime;
227  for (G4int i=0; i< myA; i++) {
228  rprime = theNucleons[i].GetPosition() -
229  factor * (theBeta*theNucleons[i].GetPosition()) * theBeta;
230  theNucleons[i].SetPosition(rprime);
231  }
232  }
233 }
234 
236 {
237  if (theBoost.e() !=0 ) {
238  G4ThreeVector beta = theBoost.vect()/theBoost.e();
239  DoLorentzContraction(beta);
240  }
241 }
242 
243 
244 
246 {
247  G4ThreeVector center;
248 
249  for (G4int i=0; i<myA; i++ )
250  {
251  center+=theNucleons[i].GetPosition();
252  }
253  center /= -myA;
254  DoTranslation(center);
255 }
256 
258 {
259  G4ThreeVector tempV;
260  for (G4int i=0; i<myA; i++ )
261  {
262  tempV = theNucleons[i].GetPosition() + theShift;
263  theNucleons[i].SetPosition(tempV);
264  }
265 }
266 
268 {
269  return theDensity;
270 }
271 
272 //----------------------- private Implementation Methods-------------
273 
275 {
276  G4int protons=0,nucleons=0;
277 
278  while (nucleons < myA ) /* Loop checking, 30-Oct-2015, G.Folger */
279  {
280  if ( protons < myZ && G4UniformRand() < (G4double)(myZ-protons)/(G4double)(myA-nucleons) )
281  {
282  protons++;
283  theNucleons[nucleons++].SetParticleType(G4Proton::Proton());
284  }
285  else if ( (nucleons-protons) < (myA-myZ) )
286  {
287  theNucleons[nucleons++].SetParticleType(G4Neutron::Neutron());
288  }
289  else G4cout << "G4Fancy3DNucleus::ChooseNucleons not efficient" << G4endl;
290  }
291  return;
292 }
293 
295 {
296  if( myA != 12) {
297 
298  G4int i=0;
299  G4ThreeVector aPos, delta;
300  G4bool freeplace;
301  const G4double nd2=sqr(nucleondistance);
302  G4double maxR=GetNuclearRadius(0.001); // there are no nucleons at a
303  // relative Density of 0.01
304  G4int jr=0;
305  G4int jx,jy;
306  G4double arand[600];
307  G4double *prand=arand;
308  places.clear(); // Reset data buffer
309  G4int interationsLeft=1000*myA;
310  while ( (i < myA) && (--interationsLeft>0)) /* Loop checking, 30-Oct-2015, G.Folger */
311  {
312  do
313  {
314  if ( jr < 3 )
315  {
316  jr=std::min(600,9*(myA - i));
317  G4RandFlat::shootArray(jr,prand);
318  //CLHEP::RandFlat::shootArray(jr, prand );
319  }
320  jx=--jr;
321  jy=--jr;
322  aPos.set((2*arand[jx]-1.), (2*arand[jy]-1.), (2*arand[--jr]-1.));
323  } while (aPos.mag2() > 1. ); /* Loop checking, 30-Oct-2015, G.Folger */
324  aPos *=maxR;
325  G4double density=theDensity->GetRelativeDensity(aPos);
326  if (G4UniformRand() < density)
327  {
328  freeplace= true;
329  std::vector<G4ThreeVector>::iterator iplace;
330  for( iplace=places.begin(); iplace!=places.end() && freeplace;++iplace)
331  {
332  delta = *iplace - aPos;
333  freeplace= delta.mag2() > nd2;
334  }
335  if ( freeplace ) {
337  // protons must at least have binding energy of CoulombBarrier, so
338  // assuming the Fermi energy corresponds to a potential, we must place these such
339  // that the Fermi Energy > CoulombBarrier
340  if (theNucleons[i].GetDefinition() == G4Proton::Proton())
341  {
342  G4double nucMass = theNucleons[i].GetDefinition()->GetPDGMass();
343  G4double eFermi= std::sqrt( sqr(pFermi) + sqr(nucMass) ) - nucMass;
344  if (eFermi <= CoulombBarrier() ) freeplace=false;
345  }
346  }
347  if ( freeplace ) {
348  theNucleons[i].SetPosition(aPos);
349  places.push_back(aPos);
350  ++i;
351  }
352  }
353  }
354  if (interationsLeft<=0) {
355  G4Exception("model/util/G4Fancy3DNucleus.cc", "mod_util001", FatalException,
356  "Problem to place nucleons");
357  }
358 
359  } else {
360  // Start insertion
361  // Alpha cluster structure of carbon nuclei, C-12, is implemented according to
362  // P. Bozek, W. Broniowski, E.R. Arriola and M. Rybczynski
363  // Phys. Rev. C90, 064902 (2014)
364  const G4double Lbase=3.05*fermi;
365  const G4double Disp=0.552; // 0.91^2*2/3 fermi^2
366  const G4double nd2=sqr(nucleondistance);
367  const G4ThreeVector Corner1=G4ThreeVector( Lbase/2., 0., 0.);
368  const G4ThreeVector Corner2=G4ThreeVector(-Lbase/2., 0., 0.);
369  const G4ThreeVector Corner3=G4ThreeVector( 0.,Lbase*0.866, 0.); // 0.866=sqrt(3)/2
370  G4ThreeVector R1;
371  R1=G4ThreeVector(G4RandGauss::shoot(0.,Disp), G4RandGauss::shoot(0.,Disp), G4RandGauss::shoot(0.,Disp))*fermi + Corner1;
372  theNucleons[0].SetPosition(R1); // First nucleon of the first He-4
373  G4int loopCounterLeft = 10000;
374  for(G4int ii=1; ii<4; ii++) // 2 - 4 nucleons of the first He-4
375  {
376  G4bool Continue;
377  do
378  {
379  R1=G4ThreeVector(G4RandGauss::shoot(0.,Disp), G4RandGauss::shoot(0.,Disp), G4RandGauss::shoot(0.,Disp))*fermi + Corner1;
380  theNucleons[ii].SetPosition(R1);
381  Continue=false;
382  for(G4int jj=0; jj < ii; jj++)
383  {
384  if( (theNucleons[ii].GetPosition() - theNucleons[jj].GetPosition()).mag2() <= nd2 ) {Continue = true; break;}
385  }
386  } while( Continue && --loopCounterLeft > 0 ); /* Loop checking, 12-Dec-2017, A.Ribon */
387  }
388  if ( loopCounterLeft <= 0 ) {
389  G4Exception("model/util/G4Fancy3DNucleus.cc", "mod_util002", FatalException,
390  "Unable to find a good position for the first alpha cluster");
391  }
392  loopCounterLeft = 10000;
393  for(G4int ii=4; ii<8; ii++) // 5 - 8 nucleons of the second He-4
394  {
395  G4bool Continue;
396  do
397  {
398  R1=G4ThreeVector(G4RandGauss::shoot(0.,Disp), G4RandGauss::shoot(0.,Disp), G4RandGauss::shoot(0.,Disp))*fermi + Corner2;
399  theNucleons[ii].SetPosition(R1);
400  Continue=false;
401  for(G4int jj=0; jj < ii; jj++)
402  {
403  if( (theNucleons[ii].GetPosition() - theNucleons[jj].GetPosition()).mag2() <= nd2 ) {Continue = true; break;}
404  }
405  } while( Continue && --loopCounterLeft > 0 ); /* Loop checking, 12-Dec-2017, A.Ribon */
406  }
407  if ( loopCounterLeft <= 0 ) {
408  G4Exception("model/util/G4Fancy3DNucleus.cc", "mod_util003", FatalException,
409  "Unable to find a good position for the second alpha cluster");
410  }
411  loopCounterLeft = 10000;
412  for(G4int ii=8; ii<12; ii++) // 9 - 12 nucleons of the third He-4
413  {
414  G4bool Continue;
415  do
416  {
417  R1=G4ThreeVector(G4RandGauss::shoot(0.,Disp), G4RandGauss::shoot(0.,Disp), G4RandGauss::shoot(0.,Disp))*fermi + Corner3;
418  theNucleons[ii].SetPosition(R1);
419  Continue=false;
420  for(G4int jj=0; jj < ii; jj++)
421  {
422  if( (theNucleons[ii].GetPosition() - theNucleons[jj].GetPosition()).mag2() <= nd2 ) {Continue = true; break;}
423  }
424  } while( Continue && --loopCounterLeft > 0 ); /* Loop checking, 12-Dec-2017, A.Ribon */
425  }
426  if ( loopCounterLeft <= 0 ) {
427  G4Exception("model/util/G4Fancy3DNucleus.cc", "mod_util004", FatalException,
428  "Unable to find a good position for the third alpha cluster");
429  }
430  G4LorentzRotation RandomRotation;
431  RandomRotation.rotateZ(2.*pi*G4UniformRand());
432  RandomRotation.rotateY(std::acos(2.*G4UniformRand()-1.));
433  // Randomly rotation of the created nucleus
435  for(G4int ii=0; ii<myA; ii++ )
436  {
437  Pos=G4LorentzVector(theNucleons[ii].GetPosition(),0.); Pos *=RandomRotation;
438  G4ThreeVector NewPos = Pos.vect();
439  theNucleons[ii].SetPosition(NewPos);
440  }
441 
442  }
443 }
444 
446 {
447  G4int i;
448  G4double density;
449 
450  // Pre-allocate buffers for filling by index
451  momentum.resize(myA, G4ThreeVector(0.,0.,0.));
452  fermiM.resize(myA, 0.*GeV);
453 
454  for (G4int ntry=0; ntry<1 ; ntry ++ )
455  {
456  for (i=0; i < myA; i++ ) // momenta for all, including last, in case we swap nucleons
457  {
458  density = theDensity->GetDensity(theNucleons[i].GetPosition());
459  fermiM[i] = theFermi.GetFermiMomentum(density);
461  if (theNucleons[i].GetDefinition() == G4Proton::Proton())
462  {
463  G4double eMax = std::sqrt(sqr(fermiM[i]) +sqr(theNucleons[i].GetDefinition()->GetPDGMass()) )
464  - CoulombBarrier();
465  if ( eMax > theNucleons[i].GetDefinition()->GetPDGMass() )
466  {
467  G4double pmax2= sqr(eMax) - sqr(theNucleons[i].GetDefinition()->GetPDGMass());
468  fermiM[i] = std::sqrt(pmax2);
469  while ( mom.mag2() > pmax2 ) /* Loop checking, 30-Oct-2015, G.Folger */
470  {
471  mom=theFermi.GetMomentum(density, fermiM[i]);
472  }
473  } else
474  {
475  //AR-21Dec2017 : emit a "JustWarning" exception instead of writing on the error stream.
476  //G4cerr << "G4Fancy3DNucleus: difficulty finding proton momentum" << G4endl;
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)",
481  "HAD_FANCY3DNUCLEUS_001", JustWarning, ed );
482  mom=G4ThreeVector(0,0,0);
483  }
484 
485  }
486  momentum[i]= mom;
487  }
488 
489  if ( ReduceSum() ) break;
490 // G4cout <<" G4FancyNucleus: iterating to find momenta: "<< ntry<< G4endl;
491  }
492 
493 // G4ThreeVector sum;
494 // for (G4int index=0; index<myA;sum+=momentum[index++])
495 // ;
496 // G4cout << "final sum / mag() " << sum << " / " << sum.mag() << G4endl;
497 
499  for ( i=0; i< myA ; i++ )
500  {
501  energy = theNucleons[i].GetParticleType()->GetPDGMass()
502  - BindingEnergy()/myA;
503  G4LorentzVector tempV(momentum[i],energy);
504  theNucleons[i].SetMomentum(tempV);
505  // GF 11-05-2011: set BindingEnergy to be T of Nucleon with p , ~ p**2/2m
506  //theNucleons[i].SetBindingEnergy(
507  // 0.5*sqr(fermiM[i])/theNucleons[i].GetParticleType()->GetPDGMass());
508  }
509 }
510 
511 
513 {
515  G4double PFermi=fermiM[myA-1];
516 
517  for (G4int i=0; i < myA-1 ; i++ )
518  { sum+=momentum[i]; }
519 
520 // check if have to do anything at all..
521  if ( sum.mag() <= PFermi )
522  {
523  momentum[myA-1]=-sum;
524  return true;
525  }
526 
527 // find all possible changes in momentum, changing only the component parallel to sum
528  G4ThreeVector testDir=sum.unit();
529  testSums.clear();
530  testSums.resize(myA-1); // Allocate block for filling below
531 
533  for (G4int aNucleon=0; aNucleon < myA-1; aNucleon++) {
534  delta = 2.*((momentum[aNucleon]*testDir)*testDir);
535 
536  testSums[aNucleon].Fill(delta, delta.mag(), aNucleon);
537  }
538 
539  std::sort(testSums.begin(), testSums.end());
540 
541 // reduce Momentum Sum until the next would be allowed.
542  G4int index=testSums.size();
543  while ( (sum-testSums[--index].Vector).mag()>PFermi && index>0) /* Loop checking, 30-Oct-2015, G.Folger */
544  {
545  // Only take one which improve, ie. don't change sign and overshoot...
546  if ( sum.mag() > (sum-testSums[index].Vector).mag() ) {
547  momentum[testSums[index].Index]-=testSums[index].Vector;
548  sum-=testSums[index].Vector;
549  }
550  }
551 
552  if ( (sum-testSums[index].Vector).mag() <= PFermi )
553  {
554  G4int best=-1;
555  G4double pBest=2*PFermi; // anything larger than PFermi
556  for ( G4int aNucleon=0; aNucleon<=index; aNucleon++)
557  {
558  // find the momentum closest to choosen momentum for last Nucleon.
559  G4double pTry=(testSums[aNucleon].Vector-sum).mag();
560  if ( pTry < PFermi
561  && std::abs(momentum[myA-1].mag() - pTry ) < pBest )
562  {
563  pBest=std::abs(momentum[myA-1].mag() - pTry );
564  best=aNucleon;
565  }
566  }
567  if ( best < 0 )
568  {
569  G4String text = "G4Fancy3DNucleus.cc: Logic error in ReduceSum()";
570  throw G4HadronicException(__FILE__, __LINE__, text);
571  }
572  momentum[testSums[best].Index]-=testSums[best].Vector;
573  momentum[myA-1]=testSums[best].Vector-sum;
574 
575  return true;
576 
577  }
578 
579  // try to compensate momentum using another Nucleon....
580  G4int swapit=-1;
581  while (swapit<myA-1) /* Loop checking, 30-Oct-2015, G.Folger */
582  {
583  if ( fermiM[++swapit] > PFermi ) break;
584  }
585  if (swapit == myA-1 ) return false;
586 
587  // Now we have a nucleon with a bigger Fermi Momentum.
588  // Exchange with last nucleon.. and iterate.
589  std::swap(theNucleons[swapit], theNucleons[myA-1]);
590  std::swap(momentum[swapit], momentum[myA-1]);
591  std::swap(fermiM[swapit], fermiM[myA-1]);
592  return ReduceSum();
593 }
594 
596 {
597  static const G4double cfactor = (1.44/1.14) * MeV;
598  return cfactor*myZ/(1.0 + G4Pow::GetInstance()->Z13(myA));
599 }