00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038 #include <algorithm>
00039
00040 #include "G4Fancy3DNucleus.hh"
00041 #include "G4Fancy3DNucleusHelper.hh"
00042 #include "G4NuclearFermiDensity.hh"
00043 #include "G4NuclearShellModelDensity.hh"
00044 #include "G4NucleiProperties.hh"
00045 #include "G4Nucleon.hh"
00046 #include "G4SystemOfUnits.hh"
00047 #include "Randomize.hh"
00048 #include "G4ios.hh"
00049 #include "G4HadronicException.hh"
00050
00051
00052 G4Fancy3DNucleus::G4Fancy3DNucleus()
00053 : myA(0), myZ(0), theNucleons(250), currentNucleon(-1), theDensity(0),
00054 nucleondistance(0.8*fermi), places(250), momentum(250), fermiM(250),
00055 testSums(250)
00056 {
00057
00058 }
00059
00060 G4Fancy3DNucleus::~G4Fancy3DNucleus()
00061 {
00062 if(theDensity) delete theDensity;
00063 }
00064
00065 #if defined(NON_INTEGER_A_Z)
00066 void G4Fancy3DNucleus::Init(G4double theA, G4double theZ)
00067 {
00068 G4int intZ = G4int(theZ);
00069 G4int intA= ( G4UniformRand()>theA-G4int(theA) ) ? G4int(theA) : G4int(theA)+1;
00070
00071 Init(intA, intZ);
00072
00073 }
00074 #endif
00075
00076 void G4Fancy3DNucleus::Init(G4int theA, G4int theZ)
00077 {
00078
00079 currentNucleon=-1;
00080 theNucleons.clear();
00081
00082 myZ = theZ;
00083 myA= theA;
00084
00085 theNucleons.resize(myA);
00086
00087
00088
00089 if(theDensity) delete theDensity;
00090 if ( myA < 17 ) {
00091 theDensity = new G4NuclearShellModelDensity(myA, myZ);
00092 } else {
00093 theDensity = new G4NuclearFermiDensity(myA, myZ);
00094 }
00095
00096 theFermi.Init(myA, myZ);
00097
00098 ChooseNucleons();
00099
00100 ChoosePositions();
00101
00102
00103
00104 ChooseFermiMomenta();
00105
00106 G4double Ebinding= BindingEnergy()/myA;
00107
00108 for (G4int aNucleon=0; aNucleon < myA; aNucleon++)
00109 {
00110 theNucleons[aNucleon].SetBindingEnergy(Ebinding);
00111 }
00112
00113
00114 return;
00115 }
00116
00117 G4bool G4Fancy3DNucleus::StartLoop()
00118 {
00119 currentNucleon=0;
00120 return (theNucleons.size()>0);
00121 }
00122
00123
00124 G4Nucleon * G4Fancy3DNucleus::GetNextNucleon()
00125 {
00126 return ( (currentNucleon>=0 && currentNucleon<myA) ?
00127 &theNucleons[currentNucleon++] : 0 );
00128 }
00129
00130 const std::vector<G4Nucleon> & G4Fancy3DNucleus::GetNucleons()
00131 {
00132 return theNucleons;
00133 }
00134
00135
00136
00137 bool G4Fancy3DNucleusHelperForSortInZ(const G4Nucleon& nuc1, const G4Nucleon& nuc2)
00138 {
00139 return nuc1.GetPosition().z() < nuc2.GetPosition().z();
00140 }
00141
00142 void G4Fancy3DNucleus::SortNucleonsIncZ()
00143 {
00144 if (theNucleons.size() < 2 ) return;
00145
00146 std::sort(theNucleons.begin(), theNucleons.end(),
00147 G4Fancy3DNucleusHelperForSortInZ);
00148 }
00149
00150 void G4Fancy3DNucleus::SortNucleonsDecZ()
00151 {
00152 if (theNucleons.size() < 2 ) return;
00153 SortNucleonsIncZ();
00154
00155 std::reverse(theNucleons.begin(), theNucleons.end());
00156 }
00157
00158
00159 G4double G4Fancy3DNucleus::BindingEnergy()
00160 {
00161 return G4NucleiProperties::GetBindingEnergy(myA,myZ);
00162 }
00163
00164
00165 G4double G4Fancy3DNucleus::GetNuclearRadius()
00166 {
00167 return GetNuclearRadius(0.5);
00168 }
00169
00170 G4double G4Fancy3DNucleus::GetNuclearRadius(const G4double maxRelativeDensity)
00171 {
00172 return theDensity->GetRadius(maxRelativeDensity);
00173 }
00174
00175 G4double G4Fancy3DNucleus::GetOuterRadius()
00176 {
00177 G4double maxradius2=0;
00178
00179 for (int i=0; i<myA; i++)
00180 {
00181 if ( theNucleons[i].GetPosition().mag2() > maxradius2 )
00182 {
00183 maxradius2=theNucleons[i].GetPosition().mag2();
00184 }
00185 }
00186 return std::sqrt(maxradius2)+nucleondistance;
00187 }
00188
00189 G4double G4Fancy3DNucleus::GetMass()
00190 {
00191 return myZ*G4Proton::Proton()->GetPDGMass() +
00192 (myA-myZ)*G4Neutron::Neutron()->GetPDGMass() -
00193 BindingEnergy();
00194 }
00195
00196
00197
00198 void G4Fancy3DNucleus::DoLorentzBoost(const G4LorentzVector & theBoost)
00199 {
00200 for (G4int i=0; i<myA; i++){
00201 theNucleons[i].Boost(theBoost);
00202 }
00203 }
00204
00205 void G4Fancy3DNucleus::DoLorentzBoost(const G4ThreeVector & theBeta)
00206 {
00207 for (G4int i=0; i<myA; i++){
00208 theNucleons[i].Boost(theBeta);
00209 }
00210 }
00211
00212 void G4Fancy3DNucleus::DoLorentzContraction(const G4ThreeVector & theBeta)
00213 {
00214 G4double factor=(1-std::sqrt(1-theBeta.mag2()))/theBeta.mag2();
00215 G4ThreeVector rprime;
00216 for (G4int i=0; i< myA; i++) {
00217 rprime = theNucleons[i].GetPosition() -
00218 factor * (theBeta*theNucleons[i].GetPosition()) * theBeta;
00219 theNucleons[i].SetPosition(rprime);
00220 }
00221 }
00222
00223 void G4Fancy3DNucleus::DoLorentzContraction(const G4LorentzVector & theBoost)
00224 {
00225 G4ThreeVector beta = theBoost.vect()/theBoost.e();
00226
00227 DoLorentzContraction(beta);
00228 }
00229
00230
00231
00232 void G4Fancy3DNucleus::CenterNucleons()
00233 {
00234 G4ThreeVector center;
00235
00236 for (G4int i=0; i<myA; i++ )
00237 {
00238 center+=theNucleons[i].GetPosition();
00239 }
00240 center /= -myA;
00241 DoTranslation(center);
00242 }
00243
00244 void G4Fancy3DNucleus::DoTranslation(const G4ThreeVector & theShift)
00245 {
00246 G4ThreeVector tempV;
00247 for (G4int i=0; i<myA; i++ )
00248 {
00249 tempV = theNucleons[i].GetPosition() + theShift;
00250 theNucleons[i].SetPosition(tempV);
00251 }
00252 }
00253
00254 const G4VNuclearDensity * G4Fancy3DNucleus::GetNuclearDensity() const
00255 {
00256 return theDensity;
00257 }
00258
00259
00260
00261 void G4Fancy3DNucleus::ChooseNucleons()
00262 {
00263 G4int protons=0,nucleons=0;
00264
00265 while (nucleons < myA )
00266 {
00267 if ( protons < myZ && G4UniformRand() < (G4double)(myZ-protons)/(G4double)(myA-nucleons) )
00268 {
00269 protons++;
00270 theNucleons[nucleons++].SetParticleType(G4Proton::Proton());
00271 }
00272 else if ( (nucleons-protons) < (myA-myZ) )
00273 {
00274 theNucleons[nucleons++].SetParticleType(G4Neutron::Neutron());
00275 }
00276 else G4cout << "G4Fancy3DNucleus::ChooseNucleons not efficient" << G4endl;
00277 }
00278 return;
00279 }
00280
00281 void G4Fancy3DNucleus::ChoosePositions()
00282 {
00283 G4int i=0;
00284 G4ThreeVector aPos, delta;
00285 G4bool freeplace;
00286 static G4double nd2 = sqr(nucleondistance);
00287 G4double maxR=GetNuclearRadius(0.001);
00288
00289 G4int jr=0;
00290 G4int jx,jy;
00291 G4double arand[600];
00292 G4double *prand=arand;
00293
00294 places.clear();
00295
00296 while ( i < myA )
00297 {
00298 do
00299 {
00300 if ( jr < 3 )
00301 {
00302 jr=std::min(600,9*(myA - i));
00303 CLHEP::RandFlat::shootArray(jr, prand );
00304 }
00305 jx=--jr;
00306 jy=--jr;
00307 aPos.set((2*arand[jx]-1.), (2*arand[jy]-1.), (2*arand[--jr]-1.));
00308 } while (aPos.mag2() > 1. );
00309 aPos *=maxR;
00310 G4double density=theDensity->GetRelativeDensity(aPos);
00311 if (G4UniformRand() < density)
00312 {
00313 freeplace= true;
00314 std::vector<G4ThreeVector>::iterator iplace;
00315 for( iplace=places.begin(); iplace!=places.end() && freeplace;++iplace)
00316 {
00317 delta = *iplace - aPos;
00318 freeplace= delta.mag2() > nd2;
00319 }
00320
00321 if ( freeplace )
00322 {
00323 G4double pFermi=theFermi.GetFermiMomentum(theDensity->GetDensity(aPos));
00324
00325
00326
00327 if (theNucleons[i].GetDefinition() == G4Proton::Proton())
00328 {
00329 G4double nucMass = theNucleons[i].GetDefinition()->GetPDGMass();
00330 G4double eFermi= std::sqrt( sqr(pFermi) + sqr(nucMass) )
00331 - nucMass;
00332 if (eFermi <= CoulombBarrier() ) freeplace=false;
00333 }
00334 }
00335 if ( freeplace )
00336 {
00337 theNucleons[i].SetPosition(aPos);
00338 places.push_back(aPos);
00339 ++i;
00340 }
00341 }
00342 }
00343
00344 }
00345
00346 void G4Fancy3DNucleus::ChooseFermiMomenta()
00347 {
00348 G4int i;
00349 G4double density;
00350
00351
00352 momentum.resize(myA, G4ThreeVector(0.,0.,0.));
00353 fermiM.resize(myA, 0.*GeV);
00354
00355 for (G4int ntry=0; ntry<1 ; ntry ++ )
00356 {
00357 for (i=0; i < myA; i++ )
00358 {
00359 density = theDensity->GetDensity(theNucleons[i].GetPosition());
00360 fermiM[i] = theFermi.GetFermiMomentum(density);
00361 G4ThreeVector mom=theFermi.GetMomentum(density);
00362 if (theNucleons[i].GetDefinition() == G4Proton::Proton())
00363 {
00364 G4double eMax = std::sqrt(sqr(fermiM[i]) +sqr(theNucleons[i].GetDefinition()->GetPDGMass()) )
00365 - CoulombBarrier();
00366 if ( eMax > theNucleons[i].GetDefinition()->GetPDGMass() )
00367 {
00368 G4double pmax2= sqr(eMax) - sqr(theNucleons[i].GetDefinition()->GetPDGMass());
00369 fermiM[i] = std::sqrt(pmax2);
00370 while ( mom.mag2() > pmax2 )
00371 {
00372 mom=theFermi.GetMomentum(density, fermiM[i]);
00373 }
00374 } else
00375 {
00376 G4cerr << "G4Fancy3DNucleus: difficulty finding proton momentum" << G4endl;
00377 mom=G4ThreeVector(0,0,0);
00378 }
00379
00380 }
00381 momentum[i]= mom;
00382 }
00383
00384 if ( ReduceSum() ) break;
00385
00386 }
00387
00388
00389
00390
00391
00392
00393 G4double energy;
00394 for ( i=0; i< myA ; i++ )
00395 {
00396 energy = theNucleons[i].GetParticleType()->GetPDGMass()
00397 - BindingEnergy()/myA;
00398 G4LorentzVector tempV(momentum[i],energy);
00399 theNucleons[i].SetMomentum(tempV);
00400
00401
00402
00403 }
00404 }
00405
00406
00407 G4bool G4Fancy3DNucleus::ReduceSum()
00408 {
00409 G4ThreeVector sum;
00410 G4double PFermi=fermiM[myA-1];
00411
00412 for (G4int i=0; i < myA-1 ; i++ )
00413 { sum+=momentum[i]; }
00414
00415
00416 if ( sum.mag() <= PFermi )
00417 {
00418 momentum[myA-1]=-sum;
00419 return true;
00420 }
00421
00422
00423 G4ThreeVector testDir=sum.unit();
00424 testSums.clear();
00425 testSums.resize(myA-1);
00426
00427 G4ThreeVector delta;
00428 for (G4int aNucleon=0; aNucleon < myA-1; aNucleon++) {
00429 delta = 2.*((momentum[aNucleon]*testDir)*testDir);
00430
00431 testSums[aNucleon].Fill(delta, delta.mag(), aNucleon);
00432 }
00433
00434 std::sort(testSums.begin(), testSums.end());
00435
00436
00437 G4int index=testSums.size();
00438 while ( (sum-testSums[--index].Vector).mag()>PFermi && index>0)
00439 {
00440
00441 if ( sum.mag() > (sum-testSums[index].Vector).mag() ) {
00442 momentum[testSums[index].Index]-=testSums[index].Vector;
00443 sum-=testSums[index].Vector;
00444 }
00445 }
00446
00447 if ( (sum-testSums[index].Vector).mag() <= PFermi )
00448 {
00449 G4int best=-1;
00450 G4double pBest=2*PFermi;
00451 for ( G4int aNucleon=0; aNucleon<=index; aNucleon++)
00452 {
00453
00454 G4double pTry=(testSums[aNucleon].Vector-sum).mag();
00455 if ( pTry < PFermi
00456 && std::abs(momentum[myA-1].mag() - pTry ) < pBest )
00457 {
00458 pBest=std::abs(momentum[myA-1].mag() - pTry );
00459 best=aNucleon;
00460 }
00461 }
00462 if ( best < 0 )
00463 {
00464 G4String text = "G4Fancy3DNucleus.cc: Logic error in ReduceSum()";
00465 throw G4HadronicException(__FILE__, __LINE__, text);
00466 }
00467 momentum[testSums[best].Index]-=testSums[best].Vector;
00468 momentum[myA-1]=testSums[best].Vector-sum;
00469
00470 return true;
00471
00472 }
00473
00474
00475 G4int swapit=-1;
00476 while (swapit<myA-1)
00477 {
00478 if ( fermiM[++swapit] > PFermi ) break;
00479 }
00480 if (swapit == myA-1 ) return false;
00481
00482
00483
00484 std::swap(theNucleons[swapit], theNucleons[myA-1]);
00485 std::swap(momentum[swapit], momentum[myA-1]);
00486 std::swap(fermiM[swapit], fermiM[myA-1]);
00487 return ReduceSum();
00488 }
00489
00490 G4double G4Fancy3DNucleus::CoulombBarrier()
00491 {
00492 G4double coulombBarrier = (1.44/1.14) * MeV * myZ / (1.0 + std::pow(G4double(myA),1./3.));
00493 return coulombBarrier;
00494 }