source: trunk/source/processes/hadronic/models/util/src/G4Fancy3DNucleus.cc@ 880

Last change on this file since 880 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 14.5 KB
Line 
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//
28// ------------------------------------------------------------
29// GEANT 4 class implementation file
30//
31// ---------------- G4Fancy3DNucleus ----------------
32// by Gunter Folger, May 1998.
33// class for a 3D nucleus, arranging nucleons in space and momentum.
34// ------------------------------------------------------------
35
36#include "G4Fancy3DNucleus.hh"
37#include "G4NuclearFermiDensity.hh"
38#include "G4NuclearShellModelDensity.hh"
39#include "G4NucleiPropertiesTable.hh"
40#include "Randomize.hh"
41#include "G4ios.hh"
42#include <algorithm>
43#include "G4HadronicException.hh"
44
45G4Fancy3DNucleus::G4Fancy3DNucleus()
46 : nucleondistance(0.8*fermi)
47{
48 theDensity=0;
49 theNucleons=0;
50 currentNucleon=-1;
51 myA=0;
52 myZ=0;
53//G4cout <<"G4Fancy3DNucleus::G4Fancy3DNucleus()"<<G4endl;
54}
55
56G4Fancy3DNucleus::~G4Fancy3DNucleus()
57{
58 if(theNucleons) delete [] theNucleons;
59 if(theDensity) delete theDensity;
60}
61
62
63void G4Fancy3DNucleus::Init(G4double theA, G4double theZ)
64{
65// G4cout << "G4Fancy3DNucleus::Init(theA, theZ) called"<<G4endl;
66 currentNucleon=-1;
67 if(theNucleons) delete [] theNucleons;
68
69// this was delected already:
70// std::for_each(theRWNucleons.begin(), theRWNucleons.end(), DeleteNucleon());
71 theRWNucleons.clear();
72
73 myZ = G4int(theZ);
74 myA= ( G4UniformRand()>theA-G4int(theA) ) ? G4int(theA) : G4int(theA)+1;
75
76 theNucleons = new G4Nucleon[myA];
77
78// G4cout << "myA, myZ" << myA << ", " << myZ << G4endl;
79
80 if(theDensity) delete theDensity;
81 if ( myA < 17 ) {
82 theDensity = new G4NuclearShellModelDensity(myA, myZ);
83 } else {
84 theDensity = new G4NuclearFermiDensity(myA, myZ);
85 }
86
87 theFermi.Init(myA, myZ);
88
89 ChooseNucleons();
90
91 ChoosePositions();
92
93// CenterNucleons(); // This would introduce a bias
94
95 ChooseFermiMomenta();
96
97 G4double Ebinding= BindingEnergy()/myA;
98
99 for (G4int aNucleon=0; aNucleon < myA; aNucleon++)
100 {
101 theNucleons[aNucleon].SetBindingEnergy(Ebinding);
102 }
103
104
105 return;
106}
107
108G4bool G4Fancy3DNucleus::StartLoop()
109{
110 currentNucleon=0;
111 return theNucleons;
112}
113
114G4Nucleon * G4Fancy3DNucleus::GetNextNucleon()
115{
116 return ( currentNucleon>=0 && currentNucleon<myA ) ?
117 theNucleons+currentNucleon++ : 0;
118}
119
120
121const std::vector<G4Nucleon *> & G4Fancy3DNucleus::GetNucleons()
122{
123 if ( theRWNucleons.size()==0 )
124 {
125 for (G4int i=0; i< myA; i++)
126 {
127 theRWNucleons.push_back(theNucleons+i);
128 }
129 }
130 return theRWNucleons;
131}
132
133G4double G4Fancy3DNucleus::BindingEnergy()
134{
135 return G4NucleiPropertiesTable::GetBindingEnergy(myZ,myA);
136}
137
138G4double G4Fancy3DNucleus::GetNuclearRadius()
139{
140 return GetNuclearRadius(0.5);
141}
142
143G4double G4Fancy3DNucleus::GetNuclearRadius(const G4double maxRelativeDensity)
144{
145 return theDensity->GetRadius(maxRelativeDensity);
146}
147
148G4double G4Fancy3DNucleus::GetOuterRadius()
149{
150 G4double maxradius2=0;
151
152 for (int i=0; i<myA; i++)
153 {
154 if ( theNucleons[i].GetPosition().mag2() > maxradius2 )
155 {
156 maxradius2=theNucleons[i].GetPosition().mag2();
157 }
158 }
159 return std::sqrt(maxradius2)+nucleondistance;
160}
161
162G4double G4Fancy3DNucleus::GetMass()
163{
164 return myZ*G4Proton::Proton()->GetPDGMass() +
165 (myA-myZ)*G4Neutron::Neutron()->GetPDGMass() -
166 BindingEnergy();
167}
168
169
170
171void G4Fancy3DNucleus::DoLorentzBoost(const G4LorentzVector & theBoost)
172{
173 for (G4int i=0; i<myA; i++){
174 theNucleons[i].Boost(theBoost);
175 }
176}
177
178void G4Fancy3DNucleus::DoLorentzBoost(const G4ThreeVector & theBeta)
179{
180 for (G4int i=0; i<myA; i++){
181 theNucleons[i].Boost(theBeta);
182 }
183}
184
185void G4Fancy3DNucleus::DoLorentzContraction(const G4ThreeVector & theBeta)
186{
187 G4double factor=(1-std::sqrt(1-theBeta.mag2()))/theBeta.mag2(); // (gamma-1)/gamma/beta**2
188 for (G4int i=0; i< myA; i++)
189 {
190 G4ThreeVector rprime=theNucleons[i].GetPosition() -
191 factor * (theBeta*theNucleons[i].GetPosition()) *
192 // theNucleons[i].GetPosition();
193 theBeta;
194 theNucleons[i].SetPosition(rprime);
195 }
196}
197
198void G4Fancy3DNucleus::DoLorentzContraction(const G4LorentzVector & theBoost)
199{
200 G4ThreeVector beta= 1/theBoost.e() * theBoost.vect();
201 // DoLorentzBoost(beta);
202 DoLorentzContraction(beta);
203}
204
205
206
207void G4Fancy3DNucleus::CenterNucleons()
208{
209 G4ThreeVector center;
210
211 for (G4int i=0; i<myA; i++ )
212 {
213 center+=theNucleons[i].GetPosition();
214 }
215 center *= -1./myA;
216 DoTranslation(center);
217}
218
219void G4Fancy3DNucleus::DoTranslation(const G4ThreeVector & theShift)
220{
221 for (G4int i=0; i<myA; i++ )
222 {
223 G4ThreeVector tempV = theNucleons[i].GetPosition() + theShift;
224 theNucleons[i].SetPosition(tempV);
225 }
226}
227
228const G4VNuclearDensity * G4Fancy3DNucleus::GetNuclearDensity() const
229{
230 return theDensity;
231}
232
233//----------------------- private Implementation Methods-------------
234
235void G4Fancy3DNucleus::ChooseNucleons()
236{
237 G4int protons=0,nucleons=0;
238
239 while (nucleons < myA )
240 {
241 if ( protons < myZ && G4UniformRand() < (G4double)(myZ-protons)/(G4double)(myA-nucleons) )
242 {
243 protons++;
244 theNucleons[nucleons++].SetParticleType(G4Proton::Proton());
245 }
246 else if ( (nucleons-protons) < (myA-myZ) )
247 {
248 theNucleons[nucleons++].SetParticleType(G4Neutron::Neutron());
249 }
250 else G4cout << "G4Fancy3DNucleus::ChooseNucleons not efficient" << G4endl;
251 }
252 return;
253}
254
255void G4Fancy3DNucleus::ChoosePositions()
256{
257 G4int i=0;
258 G4ThreeVector aPos, delta;
259 std::vector<G4ThreeVector> places;
260 places.reserve(myA);
261 G4bool freeplace;
262 static G4double nd2 = sqr(nucleondistance);
263 G4double maxR=GetNuclearRadius(0.01); // there are no nucleons at a
264 // relative Density of 0.01
265 G4int jr=0;
266 G4int jx,jy;
267 G4double arand[600];
268 G4double *prand=arand;
269// G4int Attempt=0;
270 while ( i < myA )
271 {
272 do
273 {
274// ++Attempt;
275 if ( jr < 3 )
276 {
277 jr=std::min(600,9*(myA - i));
278 CLHEP::RandFlat::shootArray(jr, prand );
279 }
280 jx=--jr;
281 jy=--jr;
282 aPos=G4ThreeVector( (2*arand[jx]-1.),
283 (2*arand[jy]-1.),
284 (2*arand[--jr]-1.));
285 } while (aPos.mag2() > 1. );
286 aPos *=maxR;
287 G4double density=theDensity->GetRelativeDensity(aPos);
288 if (G4UniformRand() < density)
289 {
290 freeplace= true;
291 std::vector<G4ThreeVector>::iterator iplace;
292 for( iplace=places.begin(); iplace!=places.end() && freeplace;++iplace)
293 {
294 delta = *iplace - aPos;
295 freeplace= delta.mag2() > nd2;
296 }
297
298 if ( freeplace )
299 {
300 G4double pFermi=theFermi.GetFermiMomentum(theDensity->GetDensity(aPos));
301 // protons must at least have binding energy of CoulombBarrier, so
302 // assuming the Fermi energy corresponds to a potential, we must place these such
303 // that the Fermi Energy > CoulombBarrier
304 if (theNucleons[i].GetDefinition() == G4Proton::Proton())
305 {
306 G4double eFermi= std::sqrt( sqr(pFermi) + sqr(theNucleons[i].GetDefinition()->GetPDGMass()) )
307 - theNucleons[i].GetDefinition()->GetPDGMass();
308 if (eFermi <= CoulombBarrier() ) freeplace=false;
309 }
310 }
311 if ( freeplace )
312 {
313 theNucleons[i].SetPosition(aPos);
314 places.push_back(aPos);
315 ++i;
316 }
317 }
318 }
319// G4cout << "Att " << myA << " " << Attempt << G4endl;
320
321}
322
323void G4Fancy3DNucleus::ChooseFermiMomenta()
324{
325 G4int i;
326 G4double density;
327 G4ThreeVector * momentum=new G4ThreeVector[myA];
328
329 G4double * fermiM=new G4double[myA];
330
331 for (G4int ntry=0; ntry<1 ; ntry ++ )
332 {
333 for (i=0; i < myA; i++ ) // momenta for all, including last, in case we swap nucleons
334 {
335 density = theDensity->GetDensity(theNucleons[i].GetPosition());
336 fermiM[i] = theFermi.GetFermiMomentum(density);
337 G4ThreeVector mom=theFermi.GetMomentum(density);
338 if (theNucleons[i].GetDefinition() == G4Proton::Proton())
339 {
340 G4double eMax = std::sqrt(sqr(fermiM[i]) +sqr(theNucleons[i].GetDefinition()->GetPDGMass()) )
341 - CoulombBarrier();
342 if ( eMax > theNucleons[i].GetDefinition()->GetPDGMass() )
343 {
344 G4double pmax2= sqr(eMax) - sqr(theNucleons[i].GetDefinition()->GetPDGMass());
345 fermiM[i] = std::sqrt(pmax2);
346 while ( mom.mag2() > pmax2 )
347 {
348 mom=theFermi.GetMomentum(density, fermiM[i]);
349 }
350 } else
351 {
352 G4cerr << "G4Fancy3DNucleus: difficulty finding proton momentum" << G4endl;
353 mom=0;
354 }
355
356 }
357 momentum[i]= mom;
358 }
359
360 if (ReduceSum(momentum,fermiM) )
361 break;
362// G4cout <<" G4FancyNucleus: iterating to find momenta: "<< ntry<< G4endl;
363 }
364
365// G4ThreeVector sum;
366// for (G4int index=0; index<myA;sum+=momentum[index++])
367// ;
368// G4cout << "final sum / mag() " << sum << " / " << sum.mag() << G4endl;
369
370
371 G4double energy;
372 for ( i=0; i< myA ; i++ )
373 {
374 energy = theNucleons[i].GetParticleType()->GetPDGMass()
375 - BindingEnergy()/myA;
376 G4LorentzVector tempV(momentum[i],energy);
377 theNucleons[i].SetMomentum(tempV);
378 }
379
380 delete [] momentum;
381 delete [] fermiM;
382}
383
384 class G4Fancy3DNucleusHelper // Helper class
385 {
386 public:
387 G4Fancy3DNucleusHelper(const G4ThreeVector &vec,const G4double size,const G4int index)
388 : Vector(vec), Size(size), anInt(index) {}
389 int operator ==(const G4Fancy3DNucleusHelper &right) const
390 {
391 return this==&right;
392 }
393 int operator < (const G4Fancy3DNucleusHelper &right) const
394 {
395 return size()<right.size();
396 }
397 const G4ThreeVector& vector() const
398 {
399 return Vector;
400 }
401 G4double size() const
402 {
403 return Size;
404 }
405 G4int index() const
406 {
407 return anInt;
408 }
409 G4Fancy3DNucleusHelper operator =(const G4Fancy3DNucleusHelper &right)
410 {
411 Vector = right.Vector;
412 Size = right.Size;
413 anInt = right.anInt;
414 return *this;
415 }
416
417 private:
418 G4Fancy3DNucleusHelper(): Vector(0), Size(0), anInt(0) {G4cout << "def ctor for MixMasch" << G4endl;}
419 G4ThreeVector Vector;
420 G4double Size;
421 G4int anInt;
422 };
423
424G4bool G4Fancy3DNucleus::ReduceSum(G4ThreeVector * momentum, G4double *pFermiM)
425{
426 G4ThreeVector sum;
427 G4double PFermi=pFermiM[myA-1];
428
429 for (G4int i=0; i < myA-1 ; i++ )
430 { sum+=momentum[i]; }
431
432// check if have to do anything at all..
433 if ( sum.mag() <= PFermi )
434 {
435 momentum[myA-1]=-sum;
436 return true;
437 }
438
439// find all possible changes in momentum, changing only the component parallel to sum
440 G4ThreeVector testDir=sum.unit();
441 std::vector<G4Fancy3DNucleusHelper> testSums; // Sorted on delta.mag()
442
443 for ( G4int aNucleon=0; aNucleon < myA-1; aNucleon++){
444 G4ThreeVector delta=2*((momentum[aNucleon]*testDir)* testDir);
445 testSums.push_back(G4Fancy3DNucleusHelper(delta,delta.mag(),aNucleon));
446 }
447 std::sort(testSums.begin(), testSums.end());
448
449// reduce Momentum Sum until the next would be allowed.
450 G4int index=testSums.size();
451 while ( (sum-testSums[--index].vector()).mag()>PFermi && index>0)
452 {
453 // Only take one which improve, ie. don't change sign and overshoot...
454 if ( sum.mag() > (sum-testSums[index].vector()).mag() ) {
455 momentum[testSums[index].index()]-=testSums[index].vector();
456 sum-=testSums[index].vector();
457 }
458 }
459
460 if ( (sum-testSums[index].vector()).mag() <= PFermi )
461 {
462 G4int best=-1;
463 G4double pBest=2*PFermi; // anything larger than PFermi
464 for ( G4int aNucleon=0; aNucleon<=index; aNucleon++)
465 {
466 // find the momentum closest to choosen momentum for last Nucleon.
467 G4double pTry=(testSums[aNucleon].vector()-sum).mag();
468 if ( pTry < PFermi
469 && std::abs(momentum[myA-1].mag() - pTry ) < pBest )
470 {
471 pBest=std::abs(momentum[myA-1].mag() - pTry );
472 best=aNucleon;
473 }
474 }
475 if ( best < 0 )
476 {
477 G4String text = "G4Fancy3DNucleus.cc: Logic error in ReduceSum()";
478 throw G4HadronicException(__FILE__, __LINE__, text);
479 }
480 momentum[testSums[best].index()]-=testSums[best].vector();
481 momentum[myA-1]=testSums[best].vector()-sum;
482
483 testSums.clear();
484 return true;
485
486 }
487 testSums.clear();
488
489 // try to compensate momentum using another Nucleon....
490
491 G4int swapit=-1;
492 while (swapit<myA-1)
493 {
494 if ( pFermiM[++swapit] > PFermi ) break;
495 }
496 if (swapit == myA-1 ) return false;
497
498 // Now we have a nucleon with a bigger Fermi Momentum.
499 // Exchange with last nucleon.. and iterate.
500// G4cout << " Nucleon to swap with : " << swapit << G4endl;
501// G4cout << " Fermi momentum test, and better.. " << PFermi << " / "
502// << theFermi.GetFermiMomentum(density) << G4endl;
503// cout << theNucleons[swapit]<< G4endl << theNucleons[myA-1] << G4endl;
504// cout << momentum[swapit] << G4endl << momentum[myA-1] << G4endl;
505 G4Nucleon swap= theNucleons[swapit];
506 G4ThreeVector mom_swap=momentum[swapit];
507 G4double pf=pFermiM[swapit];
508 theNucleons[swapit]=theNucleons[myA-1];
509 momentum[swapit]=momentum[myA-1];
510 pFermiM[swapit]=pFermiM[myA-1];
511 theNucleons[myA-1]=swap;
512 momentum[myA-1]=mom_swap;
513 pFermiM[myA-1]=pf;
514// cout << "after swap" <<G4endl<< theNucleons[swapit] << G4endl << theNucleons[myA-1] << G4endl;
515// cout << momentum[swapit] << G4endl << momentum[myA-1] << G4endl;
516 return ReduceSum(momentum,pFermiM);
517}
518
519G4double G4Fancy3DNucleus::CoulombBarrier()
520{
521 G4double coulombBarrier = (1.44/1.14) * MeV * myZ / (1.0 + std::pow(G4double(myA),1./3.));
522 return coulombBarrier;
523}
Note: See TracBrowser for help on using the repository browser.