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

Last change on this file since 1036 was 962, checked in by garnier, 17 years ago

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