Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNAChemistryManager.cc
Go to the documentation of this file.
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// Author: Mathieu Karamitros ([email protected])
28//
29// WARNING : This class is released as a prototype.
30// It might strongly evolve or even disappear in the next releases.
31//
32// History:
33// -----------
34// 10 Oct 2011 M.Karamitros created
35//
36// -------------------------------------------------------------------
37
39
40#include "G4AutoLock.hh"
44#include "G4Electron_aq.hh"
45#include "G4GeometryManager.hh"
46#include "G4H2O.hh"
48#include "G4Molecule.hh"
49#include "G4MoleculeFinder.hh"
50#include "G4MoleculeTable.hh"
51#include "G4PhysChemIO.hh"
52#include "G4Scheduler.hh"
53#include "G4StateManager.hh"
54#include "G4SystemOfUnits.hh"
55#include "G4UIcmdWithABool.hh"
59#include "G4VITTrackHolder.hh"
60#include "G4VMoleculeCounter.hh"
62#include "G4VUserPulseInfo.hh"
63
64#include <memory>
65
66G4DNAChemistryManager* G4DNAChemistryManager::fgInstance = nullptr;
67
68G4ThreadLocal G4DNAChemistryManager::ThreadLocalData*
69 G4DNAChemistryManager::fpThreadData = nullptr;
70
72
73//------------------------------------------------------------------------------
74
75G4DNAChemistryManager::ThreadLocalData::ThreadLocalData()
76{
77 fpPhysChemIO = nullptr;
78 fThreadInitialized = false;
79}
80
81//------------------------------------------------------------------------------
82
83G4DNAChemistryManager::ThreadLocalData::~ThreadLocalData()
84{
85 fpThreadData = nullptr;
86}
87
88//------------------------------------------------------------------------------
89
90void G4DNAChemistryManager::SetPhysChemIO(std::unique_ptr<G4VPhysChemIO> pPhysChemIO)
91{
92 fpThreadData->fpPhysChemIO = std::move(pPhysChemIO);
93}
94
95//------------------------------------------------------------------------------
96
97//------------------------------------------------------------------------------
98/*
99 * The chemistry manager is shared between threads
100 * It is initialized both on the master thread and on the worker threads
101 */
102//------------------------------------------------------------------------------
104 :
105 fpChemDNADirectory(new G4UIdirectory("/chem/"))
106 , fpActivateChem(new G4UIcmdWithABool("/chem/activate", this))
107 , fpRunChem(new G4UIcmdWithAnInteger("/chem/run", this))
108 , fpSkipReactionsFromChemList(new G4UIcmdWithoutParameter("/chem/skipReactionsFromChemList", this))
109 , fpScaleForNewTemperature(new G4UIcmdWithADoubleAndUnit("/chem/temperature", this))
110 , fpInitChem(new G4UIcmdWithoutParameter("/chem/init", this))
111 ,
112 fpExcitationLevel(nullptr)
113 , fpIonisationLevel(nullptr)
114 , fpUserChemistryList(nullptr)
115{
116 fpRunChem->SetParameterName("Number of runs to execute for the chemistry module"
117 "(this works when used in standalone", true, true);
118 fpRunChem->SetDefaultValue(1);
119 fpScaleForNewTemperature->SetUnitCategory("Temperature");
120}
121
122//------------------------------------------------------------------------------
123
125{
126 if (fgInstance == nullptr)
127 {
129 if (fgInstance == nullptr) // MT : double check at initialisation
130 {
131 fgInstance = new G4DNAChemistryManager();
132 }
133 lock.unlock();
134 }
135
136 // make sure thread local data is initialized for all threads
137 if (fpThreadData == nullptr)
138 {
139 fpThreadData = new ThreadLocalData();
140 }
141
142 assert(fpThreadData != nullptr);
143
144 return fgInstance;
145}
146
147//------------------------------------------------------------------------------
148
153
154//------------------------------------------------------------------------------
155
157{
158 Clear();
159 fgInstance = nullptr;
160}
161
162//------------------------------------------------------------------------------
163
165{
166 fpIonisationLevel.reset();
167 fpExcitationLevel.reset();
168
169 if (fpUserChemistryList)
170 {
171 Deregister(*fpUserChemistryList);
172 }
173
174 fpChemDNADirectory.reset();
175 fpActivateChem.reset();
176 fpRunChem.reset();
177
178 fpSkipReactionsFromChemList.reset();
179 fpInitChem.reset();
180
181 if (fpThreadData != nullptr)
182 {
183 delete fpThreadData;
184 fpThreadData = nullptr;
185 }
186
190}
191
192//------------------------------------------------------------------------------
193
195{
197
198 if (fgInstance != nullptr)
199 {
200 G4DNAChemistryManager* pDeleteMe = fgInstance;
201 fgInstance = nullptr;
202 lock.unlock();
203 delete pDeleteMe;
204 }
205 else
206 {
207 G4cerr << "G4DNAChemistryManager already deleted" << G4endl;
208 }
209 lock.unlock();
210}
211
212//------------------------------------------------------------------------------
213
215{
216 if (requestedState == G4State_Quit)
217 {
218 if (fVerbose != 0)
219 {
220 G4cout << "G4DNAChemistryManager::Notify ---> received G4State_Quit"
221 << G4endl;
222 }
223 Clear();
224 }
225 else if (requestedState == G4State_GeomClosed)
226 {
227 fGeometryClosed = true;
228 }
229 else if (requestedState == G4State_Idle)
230 {
232 }
233
234 return true;
235}
236
237//------------------------------------------------------------------------------
238
240{
241 if (pCommand == fpActivateChem.get())
242 {
244 }
245 else if (pCommand == fpRunChem.get())
246 {
247 int nbExec = value.empty() ? 1 : G4UIcommand::ConvertToInt(value);
248 for (int i = 0 ; i < nbExec ; ++i)
249 {
250 Run();
251 }
252 }
253 else if (pCommand == fpSkipReactionsFromChemList.get())
254 {
255 fSkipReactions = true;
256 }
257 else if (pCommand == fpScaleForNewTemperature.get())
258 {
259 SetGlobalTemperature(fpScaleForNewTemperature->ConvertToDimensionedDouble(value));
260 }
261 else if (pCommand == fpInitChem.get())
262 {
263 Initialize();
265 }
266}
267
268//------------------------------------------------------------------------------
269
271{
272 if (pCommand == fpActivateChem.get())
273 {
274 return G4UIcmdWithABool::ConvertToString(fActiveChemistry);
275 }
276 if (pCommand == fpScaleForNewTemperature.get())
277 {
278 return fpScaleForNewTemperature->ConvertToStringWithBestUnit(G4MolecularConfiguration::GetGlobalTemperature());
279 }
280 if (pCommand == fpSkipReactionsFromChemList.get())
281 {
282 return G4UIcmdWithABool::ConvertToString(fSkipReactions);
283 }
284
285 return "";
286}
287
288//------------------------------------------------------------------------------
289
300
301//------------------------------------------------------------------------------
303{
304 if (!fActiveChemistry)
305 {
306 return;
307 }
308
310
311 if (!fMasterInitialized)
312 {
313 G4ExceptionDescription description;
314 description << "Global components were not initialized.";
315 G4Exception("G4DNAChemistryManager::Run", "MASTER_INIT", FatalException,
316 description);
317 }
318
319 if (!fpThreadData->fThreadInitialized)
320 {
321 G4ExceptionDescription description;
322 description << "Thread local components were not initialized.";
323 G4Exception("G4DNAChemistryManager::Run", "THREAD_INIT", FatalException,
324 description);
325 }
326
329 if (fResetCounterWhenRunEnds)
330 {
332 }
333 CloseFile();
334}
335
336//------------------------------------------------------------------------------
337
339{
340 fUseInStandalone = flag;
341}
342
343//------------------------------------------------------------------------------
344
346{
347 G4Scheduler::Instance()->SetGun(pChemGun);
348}
349
350//------------------------------------------------------------------------------
351
353{
354 //===========================================================================
355 // MT MODE
356 //===========================================================================
358 {
359 //==========================================================================
360 // ON WORKER THREAD
361 //==========================================================================
363 {
364 InitializeThread(); // Will create and initialize G4Scheduler
365 return;
366 }
367 //==========================================================================
368 // ON MASTER THREAD
369 //==========================================================================
370
373 return;
374 }
375 //===========================================================================
376 // IS NOT IN MT MODE
377 //===========================================================================
378
381 // In this case: InitializeThread is called when Run() is called
382 return;
383}
384
385//------------------------------------------------------------------------------
386
388{
389 if (fMasterInitialized)
390 {
391 return;
392 }
393
394 if (fVerbose != 0)
395 {
396 G4cout << "G4DNAChemistryManager::InitializeMaster() is called" << G4endl;
397 }
398
399
400 if (fpUserChemistryList == nullptr)
401 {
402 G4ExceptionDescription description;
403 description << "No user chemistry list has been provided.";
404 G4Exception("G4DNAChemistryManager::InitializeMaster", "NO_CHEM_LIST",
405 FatalException, description);
406 }else
407 {
408 fpUserChemistryList->ConstructDissociationChannels();
409 if (!fSkipReactions)
410 {
411 fpUserChemistryList->ConstructReactionTable(G4DNAMolecularReactionTable::GetReactionTable());
412 }
413 else
414 {
416 }
417 }
418
420 // creates a concrete object of the scheduler
421 fMasterInitialized = true;
422}
423
424//------------------------------------------------------------------------------
425
427{
428 if (!fUseInStandalone || fPhysicsTableBuilt)
429 {
430 return;
431 }
432
433 if (fVerbose != 0)
434 {
435 G4cout << "G4DNAChemistryManager: Build the physics tables for "
436 "molecule definition only."
437 << G4endl;
438 }
439
440 fpUserChemistryList->BuildPhysicsTable();
441
442 if (!fGeometryClosed)
443 {
444 if (fVerbose != 0)
445 {
446 G4cout << "G4DNAChemistryManager: Close geometry" << G4endl;
447 }
448
450 pGeomManager->OpenGeometry();
451 pGeomManager->CloseGeometry(true, true);
452 fGeometryClosed = true;
453 }
454
455 fPhysicsTableBuilt = true;
456}
457
458//------------------------------------------------------------------------------
459
461{
462 if (fpThreadData->fThreadInitialized && !fForceThreadReinitialization)
463 {
464 return;
465 }
466
467 if (fpUserChemistryList == nullptr)
468 {
469 G4ExceptionDescription description;
470 description << "No user chemistry list has been provided.";
471 G4Exception("G4DNAChemistryManager::InitializeThread", "NO_CHEM_LIST",
472 FatalException, description);
473 }else
474 {
475 HandleStandaloneInitialization();// To make coverty happy
476 fpUserChemistryList->ConstructTimeStepModel(G4DNAMolecularReactionTable::GetReactionTable());
477 }
478
479 if (fVerbose != 0)
480 {
481 G4cout << "G4DNAChemistryManager::InitializeThread() is called"
482 << G4endl;
483 }
484
486
487 fpThreadData->fThreadInitialized = true;
488
490
492}
493
494//------------------------------------------------------------------------------
495
497{
498 if (fVerbose != 0)
499 {
500 G4cout << "G4DNAChemistryManager::InitializeFile() is called"
501 << G4endl;
502 }
503
504 if (fpThreadData->fpPhysChemIO)
505 {
506 fpThreadData->fpPhysChemIO->InitializeFile();
507 }
508}
509
510//------------------------------------------------------------------------------
511
513{
514 return fgInstance != nullptr ? fgInstance->IsChemistryActivated() : false;
515}
516
517//------------------------------------------------------------------------------
518
520{
521 return fActiveChemistry;
522}
523
524//------------------------------------------------------------------------------
525
527{
528 fActiveChemistry = flag;
529}
530
531//------------------------------------------------------------------------------
532
534 std::ios_base::openmode mode)
535{
536 if (fVerbose != 0)
537 {
538 G4cout << "G4DNAChemistryManager: Write chemical stage into "
539 << output.data() << G4endl;
540 }
541
542 if (!fpThreadData->fpPhysChemIO)
543 {
544 fpThreadData->fpPhysChemIO = std::make_unique<G4PhysChemIO::FormattedText>();
545 }
546
547 fpThreadData->fpPhysChemIO->WriteInto(output, mode);
548
549}
550
551//------------------------------------------------------------------------------
552
554{
555 if (fpThreadData->fpPhysChemIO)
556 {
557 fpThreadData->fpPhysChemIO->AddEmptyLineInOutputFile();
558 }
559}
560
561//------------------------------------------------------------------------------
562
564{
565 if (fpThreadData->fpPhysChemIO)
566 {
567 fpThreadData->fpPhysChemIO->CloseFile();
568 }
569}
570
571//------------------------------------------------------------------------------
572
574{
575 if (!fpExcitationLevel)
576 {
577 fpExcitationLevel = std::make_unique<G4DNAWaterExcitationStructure>();
578 }
579 return fpExcitationLevel.get();
580}
581
582//------------------------------------------------------------------------------
583
585{
586 if (!fpIonisationLevel)
587 {
588 fpIonisationLevel = std::make_unique<G4DNAWaterIonisationStructure>();
589 }
590 return fpIonisationLevel.get();
591}
592
593//------------------------------------------------------------------------------
594
596 G4int electronicLevel,
597 const G4Track* pIncomingTrack)
598{
599 if (fpThreadData->fpPhysChemIO)
600 {
601 G4double energy = -1.;
602
603 switch (modification)
604 {
606 energy = 0.;
607 break;
608 case eExcitedMolecule:
609 energy = GetExcitationLevel()->ExcitationEnergy(electronicLevel);
610 break;
611 case eIonizedMolecule:
612 energy = GetIonisationLevel()->IonisationEnergy(electronicLevel);
613 break;
614 }
615
616 fpThreadData->fpPhysChemIO->CreateWaterMolecule(modification,
617 4 - electronicLevel,
618 energy,
619 pIncomingTrack);
620 }
621
622 if (fActiveChemistry)
623 {
624 auto pH2OMolecule = new G4Molecule(G4H2O::Definition());
625
626 switch (modification)
627 {
629 pH2OMolecule->AddElectron(5, 1);
630 break;
631 case eExcitedMolecule:
632 pH2OMolecule->ExciteMolecule(4 - electronicLevel);
633 break;
634 case eIonizedMolecule:
635 pH2OMolecule->IonizeMolecule(4 - electronicLevel);
636 break;
637 }
638
639 G4double delayedTime = 0.;
640 if(pIncomingTrack->GetUserInformation() != nullptr)
641 {
642 auto pPulseInfo = dynamic_cast<G4VUserPulseInfo*>
643 (pIncomingTrack->GetUserInformation());
644 if(pPulseInfo != nullptr)
645 {
646 delayedTime = pPulseInfo->GetDelayedTime();
647 }
648 }
649
650 G4Track* pH2OTrack = pH2OMolecule->BuildTrack(picosecond + delayedTime,
651 pIncomingTrack->GetPosition());
652
653 pH2OTrack->SetParentID(pIncomingTrack->GetTrackID());
654 pH2OTrack->SetTrackStatus(fStopButAlive);
655 pH2OTrack->SetKineticEnergy(0.);
656 PushTrack(pH2OTrack);
657 }
658}
659
660//------------------------------------------------------------------------------
661// pFinalPosition is optional
663 G4ThreeVector* pFinalPosition)
664{
665 if (fpThreadData->fpPhysChemIO)
666 {
667 fpThreadData->fpPhysChemIO->CreateSolvatedElectron(pIncomingTrack,
668 pFinalPosition);
669 }
670
671 if (fActiveChemistry)
672 {
673 G4double delayedTime = 0.;
674 if(pIncomingTrack->GetUserInformation() != nullptr)
675 {
676 auto pPulseInfo = dynamic_cast<G4VUserPulseInfo*>
677 (pIncomingTrack->GetUserInformation());
678 if(pPulseInfo != nullptr)
679 {
680 delayedTime = pPulseInfo->GetDelayedTime();
681 }
682 }
683
684 PushMolecule(std::make_unique<G4Molecule>(G4Electron_aq::Definition()),
685 picosecond + delayedTime,
686 pFinalPosition != nullptr ? *pFinalPosition : pIncomingTrack->GetPosition(),
687 pIncomingTrack->GetTrackID());
688 }
689}
690
691//------------------------------------------------------------------------------
692
693void G4DNAChemistryManager::PushMolecule(std::unique_ptr<G4Molecule> pMolecule,
694 double time,
695 const G4ThreeVector& position,
696 int parentID)
697{
698 assert(fActiveChemistry
699 && "To inject chemical species, the chemistry must be activated. "
700 "Check chemistry activation before injecting species.");
701 G4Track* pTrack = pMolecule->BuildTrack(time, position);
702 pTrack->SetTrackStatus(fAlive);
703 pTrack->SetParentID(parentID);
704 pMolecule.release();
705 PushTrack(pTrack);
706}
707
708//------------------------------------------------------------------------------
709
715
716//------------------------------------------------------------------------------
717// [[deprecated]] : chemistry list should never be nullptr
719{
720 fpUserChemistryList.reset(pChemistryList);
721 fOwnChemistryList = false;
723}
724
726{
727 fpUserChemistryList.reset(&chemistryList);
728 fOwnChemistryList = false;
730}
731
732void G4DNAChemistryManager::SetChemistryList(std::unique_ptr<G4VUserChemistryList> pChemistryList)
733{
734 fpUserChemistryList = std::move(pChemistryList);
735 fOwnChemistryList = true;
737}
738
739//------------------------------------------------------------------------------
740
742{
743 if (fpUserChemistryList.get() == &chemistryList)
744 {
745 if (!fpUserChemistryList->IsPhysicsConstructor() || fOwnChemistryList)
746 {
747 fpUserChemistryList.reset();
748 }
749
750 fpUserChemistryList.release();
751 }
752}
753
754//------------------------------------------------------------------------------
755
760
761//------------------------------------------------------------------------------
762
764{
765 fVerbose = verbose;
766}
767
768//------------------------------------------------------------------------------
769
771{
772 return fResetCounterWhenRunEnds;
773}
774
775//------------------------------------------------------------------------------
776
778{
779 fResetCounterWhenRunEnds = resetCounterWhenRunEnds;
780}
781
782//------------------------------------------------------------------------------
783
785{
786 fPhysicsTableBuilt = false;
787}
788
789//------------------------------------------------------------------------------
790
792{
793 fMasterInitialized = false;
795}
796
797//------------------------------------------------------------------------------
798
800{
801 fForceThreadReinitialization = true;
802}
803
804//------------------------------------------------------------------------------
805
807{
808 fpThreadData->fThreadInitialized = false;
809}
G4ApplicationState
@ G4State_Idle
@ G4State_Quit
@ G4State_GeomClosed
G4Mutex chemManExistence
ElectronicModification
@ eIonizedMolecule
@ eDissociativeAttachment
@ eExcitedMolecule
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
std::mutex G4Mutex
@ fAlive
@ fStopButAlive
double G4double
Definition G4Types.hh:83
bool G4bool
Definition G4Types.hh:86
int G4int
Definition G4Types.hh:85
G4GLOB_DLL std::ostream G4cerr
#define G4endl
Definition G4ios.hh:67
G4GLOB_DLL std::ostream G4cout
void SetPhysChemIO(std::unique_ptr< G4VPhysChemIO > pPhysChemIO)
G4String GetCurrentValue(G4UIcommand *pCommand) override
G4bool IsCounterResetWhenRunEnds() const
void CreateSolvatedElectron(const G4Track *, G4ThreeVector *pFinalPosition=nullptr)
static G4DNAChemistryManager * GetInstanceIfExists()
static G4DNAChemistryManager * Instance()
void UseAsStandalone(G4bool flag)
void SetChemistryList(G4VUserChemistryList &)
void PushMolecule(std::unique_ptr< G4Molecule > pMolecule, G4double time, const G4ThreeVector &position, G4int parentID)
void ResetCounterWhenRunEnds(G4bool resetCounterWhenRunEnds)
void SetGlobalTemperature(G4double temperatureKelvin)
void SetGun(G4ITGun *pChemSpeciesGun)
Inject custom species to the simulation.
G4DNAWaterIonisationStructure * GetIonisationLevel()
void Deregister(G4VUserChemistryList &)
G4DNAWaterExcitationStructure * GetExcitationLevel()
void SetNewValue(G4UIcommand *, G4String) override
G4bool Notify(G4ApplicationState requestedState) override
void SetVerbose(G4int verbose)
void WriteInto(const G4String &, std::ios_base::openmode mode=std::ios_base::out)
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
static G4DNAMolecularReactionTable * GetReactionTable()
static G4DNAMolecularReactionTable * Instance()
void ScaleReactionRateForNewTemperature(double temp_K)
static G4Electron_aq * Definition()
static G4GeometryManager * GetInstance()
G4bool CloseGeometry(G4bool pOptimise=true, G4bool verbose=false, G4VPhysicalVolume *vol=nullptr)
void OpenGeometry(G4VPhysicalVolume *vol=nullptr)
static G4H2O * Definition()
Definition G4H2O.cc:42
void Push(G4Track *) override
static G4ITTrackHolder * Instance()
static void SetGlobalTemperature(G4double)
void Finalize(G4MoleculeDefinition *)
void PrepareMolecularConfiguration()
static G4MoleculeTable * Instance()
void SetGun(G4ITGun *) override
static G4Scheduler * Instance()
void Process() override
void Initialize() override
void SetTrackStatus(const G4TrackStatus aTrackStatus)
G4int GetTrackID() const
const G4ThreeVector & GetPosition() const
G4VUserTrackInformation * GetUserInformation() const
void SetKineticEnergy(const G4double aValue)
void SetParentID(const G4int aValue)
static G4bool GetNewBoolValue(const char *paramString)
static G4String ConvertToString(G4bool boolVal)
static G4int ConvertToInt(const char *st)
static void InitializeInstance()
static void DeleteInstance()
static G4VMoleculeCounter * Instance()
virtual void ResetCounter()=0
G4bool IsWorkerThread()
G4bool IsMultithreadedApplication()
G4bool IsMasterThread()
#define G4ThreadLocal
Definition tls.hh:77