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

Last change on this file since 1347 was 1340, checked in by garnier, 15 years ago

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