Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // 28 // Hadronic Process: Nuclear De-excitations 29 // by V. Lara 30 // 31 // Modified: 32 // 25.07.08 I.Pshenichnov (in collaboration wi 33 // Mishustin (FIAS, Frankfurt, INR, M 34 // Moscow, pshenich@fias.uni-frankfur 35 36 #include <numeric> 37 38 #include "G4StatMFChannel.hh" 39 #include "G4PhysicalConstants.hh" 40 #include "G4HadronicException.hh" 41 #include "Randomize.hh" 42 #include "G4Pow.hh" 43 #include "G4Exp.hh" 44 #include "G4RandomDirection.hh" 45 46 G4StatMFChannel::G4StatMFChannel() : 47 _NumOfNeutralFragments(0), 48 _NumOfChargedFragments(0) 49 { 50 Pos.resize(8); 51 Vel.resize(8); 52 Accel.resize(8); 53 } 54 55 G4StatMFChannel::~G4StatMFChannel() 56 { 57 if (!_theFragments.empty()) { 58 std::for_each(_theFragments.begin(),_theFr 59 DeleteFragment()); 60 } 61 } 62 63 G4bool G4StatMFChannel::CheckFragments(void) 64 { 65 std::deque<G4StatMFFragment*>::iterator i; 66 for (i = _theFragments.begin(); 67 i != _theFragments.end(); ++i) 68 { 69 G4int A = (*i)->GetA(); 70 G4int Z = (*i)->GetZ(); 71 if ( (A > 1 && (Z > A || Z <= 0)) || (A= 72 } 73 return true; 74 } 75 76 void G4StatMFChannel::CreateFragment(G4int A, 77 // Create a new fragment. 78 // Fragments are automatically sorted: first c 79 // then neutral ones. 80 { 81 if (Z <= 0.5) { 82 _theFragments.push_back(new G4StatMFFragme 83 _NumOfNeutralFragments++; 84 } else { 85 _theFragments.push_front(new G4StatMFFragm 86 _NumOfChargedFragments++; 87 } 88 89 return; 90 } 91 92 G4double G4StatMFChannel::GetFragmentsCoulombE 93 { 94 G4double Coulomb = 95 std::accumulate(_theFragments.begin(),_the 96 0.0, 97 [](const G4double& running 98 G4StatMFFragment*& frag 99 { 100 return running_total + f 101 } ); 102 // G4double Coulomb = 0.0; 103 // for (unsigned int i = 0;i < _theFrag 104 // Coulomb += _theFragments[i]->GetCoulom 105 return Coulomb; 106 } 107 108 G4double G4StatMFChannel::GetFragmentsEnergy(G 109 { 110 G4double Energy = 0.0; 111 112 G4double TranslationalEnergy = 1.5*T*_theFra 113 114 std::deque<G4StatMFFragment*>::const_iterato 115 for (i = _theFragments.begin(); i != _theFra 116 { 117 Energy += (*i)->GetEnergy(T); 118 } 119 return Energy + TranslationalEnergy; 120 } 121 122 G4FragmentVector * G4StatMFChannel::GetFragmen 123 G4int anZ, 124 G4double T) 125 { 126 // calculate momenta of charged fragments 127 CoulombImpulse(anA,anZ,T); 128 129 // calculate momenta of neutral fragments 130 FragmentsMomenta(_NumOfNeutralFragments, _Nu 131 132 G4FragmentVector * theResult = new G4Fragmen 133 std::deque<G4StatMFFragment*>::iterator i; 134 for (i = _theFragments.begin(); i != _theFra 135 theResult->push_back((*i)->GetFragment(T)) 136 137 return theResult; 138 } 139 140 void G4StatMFChannel::CoulombImpulse(G4int anA 141 // Aafter breakup, fragments fly away under Co 142 // This method calculates asymptotic fragments 143 { 144 // First, we have to place the fragments ins 145 PlaceFragments(anA); 146 147 // Second, we sample initial charged fragmen 148 // _NumOfChargedFragments charged fragments 149 // of the vector _theFragments (i.e. 0) 150 FragmentsMomenta(_NumOfChargedFragments, 0, 151 152 // Third, we have to figure out the asymptot 153 // For taht we have to solve equations of mo 154 SolveEqOfMotion(anA,anZ,T); 155 156 return; 157 } 158 159 void G4StatMFChannel::PlaceFragments(G4int anA 160 // This gives the position of fragments at the 161 // Fragments positions are sampled inside prol 162 { 163 G4Pow* g4calc = G4Pow::GetInstance(); 164 const G4double R0 = G4StatMFParameters::Getr 165 G4double Rsys = 2.0*R0*g4calc->Z13(anA); 166 167 G4bool TooMuchIterations; 168 do 169 { 170 TooMuchIterations = false; 171 172 // Sample the position of the first frag 173 G4double R = (Rsys - R0*g4calc->Z13(_the 174 g4calc->A13(G4UniformRand()); 175 _theFragments[0]->SetPosition(R*G4Random 176 177 178 // Sample the position of the remaining 179 G4bool ThereAreOverlaps = false; 180 std::deque<G4StatMFFragment*>::iterator 181 for (i = _theFragments.begin()+1; i != _ 182 { 183 G4int counter = 0; 184 do 185 { 186 R = (Rsys - R0*g4calc->Z13((*i)->GetA( 187 (*i)->SetPosition(R*G4RandomDirection( 188 189 // Check that there are not overlappin 190 std::deque<G4StatMFFragment*>::iterato 191 for (j = _theFragments.begin(); j != i 192 { 193 G4ThreeVector FragToFragVector = 194 (*i)->GetPosition() - (*j)->GetPositio 195 G4double Rmin = R0*(g4calc->Z13((*i)->Ge 196 g4calc->Z13((*j)->GetA())); 197 if ( (ThereAreOverlaps = (FragToFragVect 198 { break; } 199 } 200 counter++; 201 // Loop checking, 05-Aug-2015, Vladimi 202 } while (ThereAreOverlaps && counter < 1 203 204 if (counter >= 1000) 205 { 206 TooMuchIterations = true; 207 break; 208 } 209 } 210 // Loop checking, 05-Aug-2015, Vladimir 211 } while (TooMuchIterations); 212 return; 213 } 214 215 void G4StatMFChannel::FragmentsMomenta(G4int N 216 G4double T) 217 // Calculate fragments momenta at the breakup 218 // Fragment kinetic energies are calculated ac 219 // Boltzmann distribution at given temperature 220 // NF is number of fragments 221 // idx is index of first fragment 222 { 223 G4double KinE = 1.5*T*NF; 224 G4ThreeVector p(0.,0.,0.); 225 226 if (NF <= 0) return; 227 else if (NF == 1) 228 { 229 // We have only one fragment to deal wit 230 p = std::sqrt(2.0*_theFragments[idx]->Ge 231 *G4RandomDirection(); 232 _theFragments[idx]->SetMomentum(p); 233 } 234 else if (NF == 2) 235 { 236 // We have only two fragment to deal wit 237 G4double M1 = _theFragments[idx]->GetNuc 238 G4double M2 = _theFragments[idx+1]->GetN 239 p = std::sqrt(2.0*KinE*(M1*M2)/(M1+M2))* 240 _theFragments[idx]->SetMomentum(p); 241 _theFragments[idx+1]->SetMomentum(-p); 242 } 243 else 244 { 245 // We have more than two fragments 246 G4double AvailableE; 247 G4int i1,i2; 248 G4double SummedE; 249 G4ThreeVector SummedP(0.,0.,0.); 250 do 251 { 252 // Fisrt sample momenta of NF-2 fragments 253 // according to Boltzmann distribution 254 AvailableE = 0.0; 255 SummedE = 0.0; 256 SummedP.setX(0.0);SummedP.setY(0.0);Summed 257 for (G4int i = idx; i < idx+NF-2; ++i) 258 { 259 G4double E; 260 G4double RandE; 261 do 262 { 263 E = 9.0*G4UniformRand(); 264 RandE = std::sqrt(0.5/E)*G4Exp(E-0.5)*G4 265 } 266 // Loop checking, 05-Aug-2015, Vladimi 267 while (RandE > 1.0); 268 E *= T; 269 p = std::sqrt(2.0*E*_theFragments[i]-> 270 *G4RandomDirection(); 271 _theFragments[i]->SetMomentum(p); 272 SummedE += E; 273 SummedP += p; 274 } 275 // Calculate momenta of last two fragments 276 // that constraints are satisfied 277 i1 = idx+NF-2; // before last fragment in 278 i2 = idx+NF-1; // last fragment index 279 p = -SummedP; 280 AvailableE = KinE - SummedE; 281 // Available Kinetic Energy should be shar 282 } 283 // Loop checking, 05-Aug-2015, Vladimir 284 while (AvailableE <= p.mag2()/(2.0*(_the 285 _theFragments[i2]->GetNuclearMass( 286 G4double H = 1.0 + _theFragments[i2]->Ge 287 /_theFragments[i1]->GetNuclearMass(); 288 G4double CTM12 = H*(1.0 - 2.0*_theFragme 289 *AvailableE/p.mag2()); 290 G4double CosTheta1; 291 G4double Sign; 292 293 if (CTM12 > 1.) {CosTheta1 = 1.;} 294 else { 295 do 296 { 297 do 298 { 299 CosTheta1 = 1.0 - 2.0*G4UniformRand(); 300 } 301 // Loop checking, 05-Aug-2015, Vladimir 302 while (CosTheta1*CosTheta1 < CTM12); 303 } 304 // Loop checking, 05-Aug-2015, Vladimir Ivan 305 while (CTM12 >= 0.0 && CosTheta1 < 0.0); 306 } 307 308 if (CTM12 < 0.0) Sign = 1.0; 309 else if (G4UniformRand() <= 0.5) Sign = 310 else Sign = 1.0; 311 312 G4double P1 = (p.mag()*CosTheta1+Sign*st 313 *(CosTheta1*CosTheta1-CTM12) 314 G4double P2 = std::sqrt(P1*P1+p.mag2() - 315 G4double Phi = twopi*G4UniformRand(); 316 G4double SinTheta1 = std::sqrt(1.0 - Cos 317 G4double CosPhi1 = std::cos(Phi); 318 G4double SinPhi1 = std::sin(Phi); 319 G4double CosPhi2 = -CosPhi1; 320 G4double SinPhi2 = -SinPhi1; 321 G4double CosTheta2 = (p.mag2() + P2*P2 - 322 G4double SinTheta2 = 0.0; 323 if (CosTheta2 > -1.0 && CosTheta2 < 1.0) 324 SinTheta2 = std::sqrt(1.0 - CosTheta2*CosThe 325 } 326 G4ThreeVector p1(P1*SinTheta1*CosPhi1,P1 327 G4ThreeVector p2(P2*SinTheta2*CosPhi2,P2 328 G4ThreeVector b(1.0,0.0,0.0); 329 330 p1 = RotateMomentum(p,b,p1); 331 p2 = RotateMomentum(p,b,p2); 332 333 SummedP += p1 + p2; 334 SummedE += p1.mag2()/(2.0*_theFragments[ 335 p2.mag2()/(2.0*_theFragments[i2]->GetNuclear 336 337 _theFragments[i1]->SetMomentum(p1); 338 _theFragments[i2]->SetMomentum(p2); 339 340 } 341 return; 342 } 343 344 void G4StatMFChannel::SolveEqOfMotion(G4int an 345 // This method will find a solution of Newton' 346 // for fragments in the self-consistent time-d 347 { 348 G4Pow* g4calc = G4Pow::GetInstance(); 349 G4double CoulombEnergy = 0.6*CLHEP::elm_coup 350 g4calc->A13(1.0+G4StatMFParameters::GetKap 351 (G4StatMFParameters::Getr0()*g4calc->Z13(a 352 if (CoulombEnergy <= 0.0) return; 353 354 G4int Iterations = 0; 355 G4double TimeN = 0.0; 356 G4double TimeS = 0.0; 357 G4double DeltaTime = 10.0; 358 if (_NumOfChargedFragments > (G4int)Pos.size 359 Pos.resize(_NumOfChargedFragments); 360 Vel.resize(_NumOfChargedFragments); 361 Accel.resize(_NumOfChargedFragments); 362 } 363 364 G4int i; 365 for (i = 0; i < _NumOfChargedFragments; ++i) 366 { 367 Vel[i] = (1.0/(_theFragments[i]->GetNucl 368 _theFragments[i]->GetMomentum(); 369 Pos[i] = _theFragments[i]->GetPosition() 370 } 371 372 G4ThreeVector distance(0.,0.,0.); 373 G4ThreeVector force(0.,0.,0.); 374 G4ThreeVector SavedVel(0.,0.,0.); 375 do { 376 for (i = 0; i < _NumOfChargedFragments; ++ 377 { 378 force.set(0.,0.,0.); 379 for (G4int j = 0; j < _NumOfChargedFragments 380 { 381 if (i != j) 382 { 383 distance = Pos[i] - Pos[j]; 384 force += (_theFragments[i]->GetZ()*_theFra 385 (distance.mag2()*distance.mag()))*dist 386 } 387 } 388 Accel[i] = CLHEP::elm_coupling*CLHEP::fermi* 389 } 390 391 TimeN = TimeS + DeltaTime; 392 393 for ( i = 0; i < _NumOfChargedFragments; + 394 { 395 SavedVel = Vel[i]; 396 Vel[i] += Accel[i]*(TimeN-TimeS); 397 Pos[i] += (SavedVel+Vel[i])*(TimeN-TimeS)*0. 398 } 399 TimeS = TimeN; 400 401 // Loop checking, 05-Aug-2015, Vladimir Iv 402 } while (Iterations++ < 100); 403 404 // Summed fragment kinetic energy 405 G4double TotalKineticEnergy = 0.0; 406 for (i = 0; i < _NumOfChargedFragments; ++i) 407 { 408 TotalKineticEnergy += _theFragments[i]-> 409 0.5*Vel[i].mag2(); 410 } 411 // Scaling of fragment velocities 412 G4double KineticEnergy = 1.5*_NumOfChargedFr 413 G4double Eta = std::pow(( CoulombEnergy + Ki 414 415 // Finally calculate fragments momenta 416 for (i = 0; i < _NumOfChargedFragments; ++i) 417 { 418 _theFragments[i]->SetMomentum((_theFragm 419 } 420 return; 421 } 422 423 G4ThreeVector G4StatMFChannel::RotateMomentum( 424 G4ThreeVector V, G4ThreeVector 425 // Rotates a 3-vector P to close momentum 426 { 427 G4ThreeVector U = Pa.unit(); 428 429 G4double Alpha1 = U * V; 430 431 G4double Alpha2 = std::sqrt(V.mag2() - Alpha 432 433 G4ThreeVector N = (1./Alpha2)*U.cross(V); 434 435 G4ThreeVector RotatedMomentum( 436 ( (V.x() - Alpha1*U.x())/Alpha2 ) * P. 437 ( (V.y() - Alpha1*U.y())/Alpha2 ) * P. 438 ( (V.z() - Alpha1*U.z())/Alpha2 ) * P. 439 ); 440 return RotatedMomentum; 441 } 442 443