Geant4 Cross Reference |
1 // ******************************************************************** 2 // * License and Disclaimer * 3 // * * 4 // * The Geant4 software is copyright of the Copyright Holders of * 5 // * the Geant4 Collaboration. It is provided under the terms and * 6 // * conditions of the Geant4 Software License, included in the file * 7 // * LICENSE and available at http://cern.ch/geant4/license . These * 8 // * include a list of copyright holders. * 9 // * * 10 // * Neither the authors of this software system, nor their employing * 11 // * institutes,nor the agencies providing financial support for this * 12 // * work make any representation or warranty, express or implied, * 13 // * regarding this software system or assume any liability for its * 14 // * use. Please see the license in the file LICENSE and URL above * 15 // * for the full disclaimer and the limitation of liability. * 16 // * * 17 // * This code implementation is the result of the scientific and * 18 // * technical work of the GEANT4 collaboration. * 19 // * By using, copying, modifying or distributing the software (or * 20 // * any work based on the software) you agree to acknowledge its * 21 // * use in resulting scientific publications, and indicate your * 22 // * acceptance of all terms of the Geant4 Software license. * 23 // ******************************************************************** 24 // 25 // 26 // --------------------------------------------------------------------- 27 // GEANT 4 class header file 28 // 29 // History: first implementation, based on G4DynamicParticle 30 // New dependency : G4VUserTrackInformation 31 // 32 // ---------------- G4Molecule ---------------- 33 // first design&implementation by Alfonso Mantero, 7 Apr 2009 34 // New developments Alfonso Mantero & Mathieu Karamitros 35 // Oct/Nov 2009 Class Name changed to G4Molecule 36 // Removed dependency from G4DynamicParticle 37 // New constructors : 38 // copy constructor 39 // direct ionized/excited molecule 40 // New methods : 41 // Get : name,atoms' number,nb electrons,decayChannel 42 // PrintState //To get the electronic level and the 43 // corresponding name of the excitation 44 // Kinematic : 45 // BuildTrack,GetKineticEnergy,GetDiffusionVelocity 46 // Change the way dynCharge and eNb is calculated 47 // --------------------------------------------------------------------- 48 49 #include "G4Molecule.hh" 50 #include "G4MolecularConfiguration.hh" 51 #include "Randomize.hh" 52 #include "G4PhysicalConstants.hh" 53 #include "G4SystemOfUnits.hh" 54 #include "G4Track.hh" 55 #include "G4VMoleculeCounter.hh" 56 57 using namespace std; 58 59 G4Allocator<G4Molecule>*& aMoleculeAllocator() 60 { 61 G4ThreadLocalStatic G4Allocator<G4Molecule>* _instance = nullptr; 62 return _instance; 63 } 64 65 //______________________________________________________________________________ 66 67 template<> 68 G4KDNode<G4Molecule>::~G4KDNode() { 69 fPoint->SetNode(nullptr); 70 } 71 72 //______________________________________________________________________________ 73 74 G4Molecule* GetMolecule(const G4Track& track) 75 { 76 return (G4Molecule*)(GetIT(track)); 77 } 78 79 //______________________________________________________________________________ 80 81 G4Molecule* GetMolecule(const G4Track* track) 82 { 83 return (G4Molecule*)(GetIT(track)); 84 } 85 86 //______________________________________________________________________________ 87 88 G4Molecule* G4Molecule::GetMolecule(const G4Track* track) 89 { 90 return (G4Molecule*)(GetIT(track)); 91 } 92 93 //______________________________________________________________________________ 94 95 void G4Molecule::Print() const 96 { 97 G4cout << "The user track information is a molecule" << G4endl; 98 } 99 100 //______________________________________________________________________________ 101 102 G4Molecule::G4Molecule(const G4Molecule& right) 103 : G4VUserTrackInformation("G4Molecule") 104 , G4IT(right) 105 { 106 fpMolecularConfiguration = right.fpMolecularConfiguration; 107 } 108 109 //______________________________________________________________________________ 110 111 G4Molecule& G4Molecule::operator=(const G4Molecule& right) 112 { 113 if (&right == this) return *this; 114 fpMolecularConfiguration = right.fpMolecularConfiguration; 115 return *this; 116 } 117 118 //______________________________________________________________________________ 119 120 G4bool G4Molecule::operator==(const G4Molecule& right) const 121 { 122 return fpMolecularConfiguration == right.fpMolecularConfiguration; 123 } 124 125 //______________________________________________________________________________ 126 127 G4bool G4Molecule::operator!=(const G4Molecule& right) const 128 { 129 return !(*this == right); 130 } 131 132 //______________________________________________________________________________ 133 /** The two methods below are the most called of the simulation : 134 * compare molecules in the MoleculeStackManager or in 135 * the InteractionTable 136 */ 137 138 G4bool G4Molecule::operator<(const G4Molecule& right) const 139 { 140 return fpMolecularConfiguration < right.fpMolecularConfiguration; 141 } 142 143 //______________________________________________________________________________ 144 145 G4Molecule::G4Molecule() 146 : G4VUserTrackInformation("G4Molecule") 147 { 148 fpMolecularConfiguration = nullptr; 149 } 150 151 //______________________________________________________________________________ 152 153 G4Molecule::~G4Molecule() 154 { 155 if (fpTrack != nullptr) 156 { 157 if (G4VMoleculeCounter::Instance()->InUse()) 158 { 159 G4VMoleculeCounter::Instance()-> 160 RemoveAMoleculeAtTime(fpMolecularConfiguration, 161 fpTrack->GetGlobalTime(), 162 &(fpTrack->GetPosition())); 163 } 164 fpTrack = nullptr; 165 } 166 fpMolecularConfiguration = nullptr; 167 } 168 169 //______________________________________________________________________________ 170 /** Build a molecule at ground state according to a given 171 * G4MoleculeDefinition that can be obtained from G4GenericMoleculeManager 172 */ 173 G4Molecule::G4Molecule(G4MoleculeDefinition* pMoleculeDefinition) 174 : G4VUserTrackInformation("G4Molecule") 175 { 176 fpMolecularConfiguration = G4MolecularConfiguration::GetOrCreateMolecularConfiguration(pMoleculeDefinition); 177 } 178 179 //______________________________________________________________________________ 180 181 G4Molecule::G4Molecule(G4MoleculeDefinition* pMoleculeDefinition, int charge) 182 { 183 fpMolecularConfiguration = G4MolecularConfiguration::GetOrCreateMolecularConfiguration(pMoleculeDefinition, 184 charge); 185 } 186 187 //______________________________________________________________________________ 188 /** Build a molecule at a specific excitation/ionisation state according 189 * to a ground state that can be obtained from G4GenericMoleculeManager. 190 * Put 0 in the second option if this is a ionisation. 191 */ 192 G4Molecule::G4Molecule(G4MoleculeDefinition* pMoleculeDefinition, 193 G4int OrbitalToFree, 194 G4int OrbitalToFill) 195 : G4VUserTrackInformation("G4Molecule") 196 { 197 if (pMoleculeDefinition->GetGroundStateElectronOccupancy() != nullptr) 198 { 199 G4ElectronOccupancy dynElectronOccupancy(*pMoleculeDefinition->GetGroundStateElectronOccupancy()); 200 201 if (OrbitalToFill != 0) 202 { 203 dynElectronOccupancy.RemoveElectron(OrbitalToFree - 1, 1); 204 dynElectronOccupancy.AddElectron(OrbitalToFill - 1, 1); 205 // dynElectronOccupancy.DumpInfo(); // DEBUG 206 } 207 208 if (OrbitalToFill == 0) 209 { 210 dynElectronOccupancy.RemoveElectron(OrbitalToFree - 1, 1); 211 // dynElectronOccupancy.DumpInfo(); // DEBUG 212 } 213 214 fpMolecularConfiguration = 215 G4MolecularConfiguration::GetOrCreateMolecularConfiguration( 216 pMoleculeDefinition, dynElectronOccupancy); 217 } 218 else 219 { 220 fpMolecularConfiguration = nullptr; 221 G4Exception( 222 "G4Molecule::G4Molecule(G4MoleculeDefinition* pMoleculeDefinition, " 223 "G4int OrbitalToFree, G4int OrbitalToFill)", 224 "G4Molecule_wrong_usage_of_constructor", 225 FatalErrorInArgument, 226 "If you want to use this constructor, the molecule definition has to be " 227 "first defined with electron occupancies"); 228 } 229 } 230 231 //______________________________________________________________________________ 232 /** Specific builder for water molecules to be used in Geant4-DNA, 233 * the last option Excitation is true if the molecule is excited, is 234 * false is the molecule is ionized. 235 */ 236 G4Molecule::G4Molecule(G4MoleculeDefinition* pMoleculeDefinition, 237 G4int level, 238 G4bool excitation) 239 : G4VUserTrackInformation("G4Molecule") 240 { 241 if (pMoleculeDefinition->GetGroundStateElectronOccupancy() != nullptr) 242 { 243 G4ElectronOccupancy dynElectronOccupancy(*pMoleculeDefinition->GetGroundStateElectronOccupancy()); 244 245 if (excitation) 246 { 247 dynElectronOccupancy.RemoveElectron(level, 1); 248 dynElectronOccupancy.AddElectron(5, 1); 249 // dynElectronOccupancy.DumpInfo(); // DEBUG 250 } 251 else 252 { 253 dynElectronOccupancy.RemoveElectron(level, 1); 254 // dynElectronOccupancy.DumpInfo(); // DEBUG 255 } 256 257 fpMolecularConfiguration = G4MolecularConfiguration::GetOrCreateMolecularConfiguration(pMoleculeDefinition, 258 dynElectronOccupancy); 259 } 260 else 261 { 262 fpMolecularConfiguration = nullptr; 263 G4Exception( 264 "G4Molecule::G4Molecule(G4MoleculeDefinition* pMoleculeDefinition, " 265 "G4int OrbitalToFree, G4int OrbitalToFill)", 266 "G4Molecule_wrong_usage_of_constructor", 267 FatalErrorInArgument, 268 "If you want to use this constructor, the molecule definition has to be " 269 "first defined with electron occupancies"); 270 } 271 } 272 273 //______________________________________________________________________________ 274 275 G4Molecule::G4Molecule(const G4MolecularConfiguration* pMolecularConfiguration) 276 { 277 fpMolecularConfiguration = pMolecularConfiguration; 278 } 279 280 //______________________________________________________________________________ 281 282 void G4Molecule::SetElectronOccupancy(const G4ElectronOccupancy* pElectronOcc) 283 { 284 fpMolecularConfiguration = 285 G4MolecularConfiguration::GetOrCreateMolecularConfiguration(fpMolecularConfiguration->GetDefinition(), 286 *pElectronOcc); 287 } 288 289 //______________________________________________________________________________ 290 291 void G4Molecule::ExciteMolecule(G4int excitationLevel) 292 { 293 fpMolecularConfiguration = fpMolecularConfiguration->ExciteMolecule(excitationLevel); 294 } 295 296 //______________________________________________________________________________ 297 298 void G4Molecule::IonizeMolecule(G4int ionizationLevel) 299 { 300 fpMolecularConfiguration = fpMolecularConfiguration->IonizeMolecule(ionizationLevel); 301 } 302 303 //______________________________________________________________________________ 304 305 void G4Molecule::AddElectron(G4int orbit, G4int number) 306 { 307 fpMolecularConfiguration = fpMolecularConfiguration->AddElectron(orbit, number); 308 } 309 310 //______________________________________________________________________________ 311 312 void G4Molecule::RemoveElectron(G4int orbit, G4int number) 313 { 314 fpMolecularConfiguration = 315 fpMolecularConfiguration->RemoveElectron(orbit, number); 316 } 317 318 //______________________________________________________________________________ 319 320 void G4Molecule::MoveOneElectron(G4int orbitToFree, G4int orbitToFill) 321 { 322 fpMolecularConfiguration = 323 fpMolecularConfiguration->MoveOneElectron(orbitToFree, orbitToFill); 324 } 325 326 //______________________________________________________________________________ 327 328 const G4String& G4Molecule::GetName() const 329 { 330 return fpMolecularConfiguration->GetName(); 331 } 332 333 //______________________________________________________________________________ 334 335 const G4String& G4Molecule::GetFormatedName() const 336 { 337 return fpMolecularConfiguration->GetFormatedName(); 338 } 339 340 //______________________________________________________________________________ 341 342 G4int G4Molecule::GetAtomsNumber() const 343 { 344 return fpMolecularConfiguration->GetAtomsNumber(); 345 } 346 347 //______________________________________________________________________________ 348 349 G4double G4Molecule::GetNbElectrons() const 350 { 351 return fpMolecularConfiguration->GetNbElectrons(); 352 } 353 354 //______________________________________________________________________________ 355 356 void G4Molecule::PrintState() const 357 { 358 fpMolecularConfiguration->PrintState(); 359 } 360 361 //______________________________________________________________________________ 362 363 G4Track* G4Molecule::BuildTrack(G4double globalTime, 364 const G4ThreeVector& position) 365 { 366 if (fpTrack != nullptr) 367 { 368 G4Exception("G4Molecule::BuildTrack", "Molecule001", FatalErrorInArgument, 369 "A track was already assigned to this molecule"); 370 } 371 372 // Kinetic Values 373 // Set a random direction to the molecule 374 G4double costheta = (2 * G4UniformRand() - 1); 375 G4double theta = acos(costheta); 376 G4double phi = 2 * pi * G4UniformRand(); 377 378 G4double xMomentum = cos(phi) * sin(theta); 379 G4double yMomentum = sin(theta) * sin(phi); 380 G4double zMomentum = costheta; 381 382 G4ThreeVector MomentumDirection(xMomentum, yMomentum, zMomentum); 383 G4double KineticEnergy = GetKineticEnergy(); 384 385 auto dynamicParticle = new G4DynamicParticle( 386 fpMolecularConfiguration->GetDefinition(), MomentumDirection, 387 KineticEnergy); 388 389 if (G4VMoleculeCounter::Instance()->InUse()) 390 { 391 G4VMoleculeCounter::Instance()-> 392 AddAMoleculeAtTime(fpMolecularConfiguration, 393 globalTime, 394 &(fpTrack->GetPosition())); 395 } 396 397 //Set the Track 398 fpTrack = new G4Track(dynamicParticle, globalTime, position); 399 fpTrack->SetUserInformation(this); 400 401 return fpTrack; 402 } 403 404 //______________________________________________________________________________ 405 406 G4double G4Molecule::GetKineticEnergy() const 407 { 408 //// 409 // Ideal Gaz case 410 double v = GetDiffusionVelocity(); 411 double E = (fpMolecularConfiguration->GetMass() / (c_squared)) * (v * v) / 2.; 412 //// 413 return E; 414 } 415 416 //______________________________________________________________________________ 417 418 G4double G4Molecule::GetDiffusionVelocity() const 419 { 420 double moleculeMass = fpMolecularConfiguration->GetMass() / (c_squared); 421 422 //// 423 // Different possibilities 424 //// 425 // Ideal Gaz case : Maxwell Boltzmann Distribution 426 // double sigma = k_Boltzmann * fgTemperature / mass; 427 // return G4RandGauss::shoot( 0, sigma ); 428 //// 429 // Ideal Gaz case : mean velocity from equipartition theorem 430 return sqrt(3 * k_Boltzmann * 431 G4MolecularConfiguration::GetGlobalTemperature() / moleculeMass); 432 //// 433 // Using this approximation for liquid is wrong 434 // However the brownian process avoid taking 435 // care of energy consideration and plays only 436 // with positions 437 } 438 439 //______________________________________________________________________________ 440 441 // added - to be transformed in a "Decay method" 442 const vector<const G4MolecularDissociationChannel*>* 443 G4Molecule::GetDissociationChannels() const 444 { 445 return fpMolecularConfiguration->GetDissociationChannels(); 446 } 447 448 //______________________________________________________________________________ 449 450 G4int G4Molecule::GetFakeParticleID() const 451 { 452 return fpMolecularConfiguration->GetFakeParticleID(); 453 } 454 455 //______________________________________________________________________________ 456 457 G4int G4Molecule::GetMoleculeID() const 458 { 459 return fpMolecularConfiguration->GetMoleculeID(); 460 } 461 462 //______________________________________________________________________________ 463 464 G4double G4Molecule::GetDecayTime() const 465 { 466 return fpMolecularConfiguration->GetDecayTime(); 467 } 468 469 //______________________________________________________________________________ 470 471 G4double G4Molecule::GetVanDerVaalsRadius() const 472 { 473 return fpMolecularConfiguration->GetVanDerVaalsRadius(); 474 } 475 476 //______________________________________________________________________________ 477 478 G4int G4Molecule::GetCharge() const 479 { 480 return fpMolecularConfiguration->GetCharge(); 481 } 482 483 //______________________________________________________________________________ 484 485 G4double G4Molecule::GetMass() const 486 { 487 return fpMolecularConfiguration->GetMass(); 488 } 489 490 //______________________________________________________________________________ 491 492 const G4ElectronOccupancy* G4Molecule::GetElectronOccupancy() const 493 { 494 return fpMolecularConfiguration->GetElectronOccupancy(); 495 } 496 497 //______________________________________________________________________________ 498 499 const G4MoleculeDefinition* G4Molecule::GetDefinition() const 500 { 501 return fpMolecularConfiguration->GetDefinition(); 502 } 503 504 //______________________________________________________________________________ 505 506 G4double G4Molecule::GetDiffusionCoefficient() const 507 { 508 return fpMolecularConfiguration->GetDiffusionCoefficient(); 509 } 510 511 //______________________________________________________________________________ 512 513 G4double G4Molecule::GetDiffusionCoefficient(const G4Material* pMaterial, 514 double temperature) const 515 { 516 return fpMolecularConfiguration->GetDiffusionCoefficient(pMaterial, 517 temperature); 518 } 519 520 //______________________________________________________________________________ 521 522 const G4MolecularConfiguration* G4Molecule::GetMolecularConfiguration() const 523 { 524 return fpMolecularConfiguration; 525 } 526 527 //______________________________________________________________________________ 528 529 const G4String& G4Molecule::GetLabel() const 530 { 531 return fpMolecularConfiguration->GetLabel(); 532 } 533 534 //______________________________________________________________________________ 535 536 void G4Molecule::ChangeConfigurationToLabel(const G4String& label) 537 { 538 // TODO check fpMolecularConfiguration already exists 539 // and new one as well 540 // TODO notify for stack change 541 fpMolecularConfiguration = G4MolecularConfiguration::GetMolecularConfiguration( 542 fpMolecularConfiguration->GetDefinition(), label); 543 544 assert(fpMolecularConfiguration != nullptr); 545 } 546