Geant4  10.01
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 // $Id: G4DNAChemistryManager.cc 87375 2014-12-02 08:17:28Z gcosmo $
27 //
28 // Author: Mathieu Karamitros (kara@cenbg.in2p3.fr)
29 //
30 // WARNING : This class is released as a prototype.
31 // It might strongly evolve or even disapear in the next releases.
32 //
33 // History:
34 // -----------
35 // 10 Oct 2011 M.Karamitros created
36 //
37 // -------------------------------------------------------------------
38 
39 #include "G4DNAChemistryManager.hh"
40 
41 #include "G4Scheduler.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4Molecule.hh"
44 #include "G4VITTrackHolder.hh"
45 #include "G4H2O.hh"
49 #include "G4Electron_aq.hh"
51 #include "G4MoleculeCounter.hh"
52 #include "G4VUserChemistryList.hh"
53 #include "G4AutoLock.hh"
54 #include "G4UIcmdWithABool.hh"
56 #include "G4GeometryManager.hh"
57 #include "G4StateManager.hh"
58 #include "G4MoleculeFinder.hh"
59 
60 using namespace std;
61 
66 //bool G4DNAChemistryManager::fActiveChemistry = false;
67 
70 {
71 //==============================================================================
72 /* M.K: 24/11/2014
73  * To work properly, the chemistry manager should be created and initialized on
74  * the master thread only. If the static flag fActiveChemistry is on but the
75  * chemistry manager singleton
76  */
77 //==============================================================================
78 
79 // if (/*fActiveChemistry &&*/ G4Threading::IsWorkerThread()
80 // && G4Threading::IsMultithreadedApplication())
81 // {
82 // G4Exception("G4DNAChemistryManager::G4DNAChemistryManager",
83 // "G4DNAChemistryManager_MASTER_CREATION", FatalException,
84 // "The chemistry manager should be created and initialized on the "
85 // "master thread only");
86 // }
87 
90  fWriteFile = false;
92  fMasterInitialized = false;
93  fpChemDNADirectory = new G4UIdirectory("/process/em/dna/chem/");
94  fpActivateChem = new G4UIcmdWithABool("/process/em/dna/chem/activate", this);
95  fpRunChem = new G4UIcmdWithoutParameter("/process/em/dna/chem/run", this);
96  fBuildPhysicsTable = false;
97  fGeometryClosed = false;
98  fPhysicsTableBuilt = false;
100  fFileInitialized = false;
101  fVerbose = 0;
102  fActiveChemistry = false;
103 }
104 
107 {
108  if (fgInstance == 0)
109  {
111  if (fgInstance == 0) // MT : double check at initialisation
112  {
114  }
115  lock.unlock();
116  }
117  return fgInstance;
118 }
119 
122 {
123  return fgInstance;
124 }
125 
127 {
128 // G4cout << "Deleting G4DNAChemistryManager" << G4endl;
129  Clear();
130  fgInstance = 0;
131  /*
132  * DEBUG : check that the chemistry manager has well been deregistered
133  * assert(G4StateManager::GetStateManager()->
134  * DeregisterDependent(this) == true);
135  */
136 }
137 
139 {
140  if (fpIonisationLevel)
141  {
142  delete fpIonisationLevel;
143  fpIonisationLevel = 0;
144 
145  }
146  if (fpExcitationLevel)
147  {
148  delete fpExcitationLevel;
149  fpExcitationLevel = 0;
150  }
152  {
154  {
155  delete fpUserChemistryList;
156  }
157 // else
158 // {
159 // G4cout << "G4DNAChemistryManager will not delete the chemistry list "
160 // "since it inherits from G4VPhysicsConstructor and it is then "
161 // "expected to be the responsability to the G4VModularPhysics to handle"
162 // " the chemistry list." << G4endl;
163 // }
165  }
166 
167  if (fpChemDNADirectory)
168  {
169  delete fpChemDNADirectory;
170  fpChemDNADirectory = 0;
171  }
172  if (fpActivateChem)
173  {
174  delete fpActivateChem;
175  fpActivateChem = 0;
176  }
177 
178  if(fpRunChem)
179  {
180  delete fpRunChem;
181  fpRunChem = 0;
182  }
183 
185  //G4MoleculeHandleManager::DeleteInstance();
188 }
189 
191 {
192  //G4cout << "G4DNAChemistryManager::DeleteInstance" << G4endl;
193 
195 
196  if(fgInstance)
197  {
198  G4DNAChemistryManager* deleteMe = fgInstance;
199  fgInstance = 0;
200  lock.unlock();
201  delete deleteMe;
202  }
203  else
204  {
205  G4cout << "G4DNAChemistryManager already deleted" << G4endl;
206  }
207  lock.unlock();
208 }
209 
211 {
212  if (requestedState == G4State_Quit)
213  {
214  if(fVerbose)
215  G4cout << "G4DNAChemistryManager::Notify ---> received G4State_Quit"
216  << G4endl;
217  //DeleteInstance();
218  Clear();
219  }
220 
221  else if(requestedState == G4State_GeomClosed)
222  {
223  fGeometryClosed = true;
224  }
225 
226  return true;
227 }
228 
230 {
231  if (command == fpActivateChem)
232  {
234  }
235  else if (command == fpRunChem)
236  {
237  Run();
238  }
239 }
240 
242 {
243  if (fActiveChemistry)
244  {
246 
247  if (fMasterInitialized == false)
248  {
249  G4ExceptionDescription description;
250  description << "Global components were not initialized.";
251  G4Exception("G4DNAChemistryManager::Run", "MASTER_INIT", FatalException,
252  description);
253  }
254 
255  if (fpgThreadInitialized_tl == 0)
256  {
257  G4ExceptionDescription description;
258  description << "Thread local components were not initialized.";
259  G4Exception("G4DNAChemistryManager::Run", "THREAD_INIT", FatalException,
260  description);
261  }
262 
264  CloseFile();
265  }
266 }
267 
268 void G4DNAChemistryManager::Gun(G4ITGun* gun, bool physicsTableToBuild)
269 {
270  fBuildPhysicsTable = physicsTableToBuild;
272 }
273 
275 {
276  //===========================================================================
277  // MT MODE
278  //===========================================================================
280  {
281  //==========================================================================
282  // ON WORKER THREAD
283  //==========================================================================
285  {
286  InitializeThread(); // Will create and initialize G4ITScheduler
287  return;
288  }
289  //==========================================================================
290  // ON MASTER THREAD
291  //==========================================================================
292  else
293  {
295  return;
296  }
297  }
298  //===========================================================================
299  // IS NOT IN MT MODE
300  //===========================================================================
301  else
302  {
304  // In this case: InitializeThread is called when Run() is called
305  return;
306  }
307 
308 }
309 
311 {
312  if (fMasterInitialized == false)
313  {
314  if(fVerbose)
315  {
316  G4cout << "G4DNAChemistryManager::InitializeMaster() is called" << G4endl;
317  }
318 
320  // creates a concrete object of the scheduler
321  // and track container
322 
324  {
328  fMasterInitialized = true;
329  }
330  else
331  {
332  if (fActiveChemistry)
333  {
334  G4ExceptionDescription description;
335  description << "No user chemistry list has been provided.";
336  G4Exception("G4DNAChemistryManager::InitializeMaster", "NO_CHEM_LIST",
337  FatalException, description);
338  }
339  }
340  }
341 }
342 
344 {
346  {
348  {
349  if(fVerbose)
350  {
351  G4cout << "G4DNAChemistryManager::InitializeThread() is called"
352  << G4endl;
353  }
354 
355  if (fBuildPhysicsTable && fPhysicsTableBuilt == false)
356  {
357  if(fVerbose)
358  {
359  G4cout << "G4DNAChemistryManager: Build the physics tables for "
360  "molecules."
361  << G4endl;
362  }
363 
365  if (fGeometryClosed == false)
366  {
367  if(fVerbose)
368  {
369  G4cout << "G4DNAChemistryManager: Close geometry"
370  << G4endl;
371  }
372 
374  // G4cout << "Start closing geometry." << G4endl;
375  geomManager->OpenGeometry();
376  geomManager->CloseGeometry(true, true);
377  fGeometryClosed = true;
378  }
379 
380  fPhysicsTableBuilt = true;
381  }
385 
386  fpgThreadInitialized_tl = new G4bool(true);
387  }
388  else
389  {
390  G4ExceptionDescription description;
391  description << "No user chemistry list has been provided.";
392  G4Exception("G4DNAChemistryManager::InitializeThread", "NO_CHEM_LIST",
393  FatalException, description);
394  }
395 
397  }
398 
399  InitializeFile();
400 }
401 
403 {
404  if (fpgOutput_tl == 0 || fWriteFile == false || fFileInitialized)
405  {
406  return;
407  }
408 
409  if(fVerbose)
410  {
411  G4cout << "G4DNAChemistryManager::InitializeFile() is called"
412  << G4endl;
413  }
414 
415  *fpgOutput_tl << std::setprecision(6) << std::scientific;
416  *fpgOutput_tl << setw(11) << left << "#Parent ID" << setw(10) << "Molecule"
417  << setw(14) << "Elec Modif" << setw(13) << "Energy (eV)"
418  << setw(22) << "X pos of parent [nm]" << setw(22)
419  << "Y pos of parent [nm]" << setw(22) << "Z pos of parent [nm]"
420  << setw(14) << "X pos [nm]" << setw(14) << "Y pos [nm]"
421  << setw(14) << "Z pos [nm]" << G4endl<< setw(21) << "#"
422  << setw(13) << "1)io/ex=0/1"
423  << G4endl
424  << setw(21) << "#"
425  << setw(13) << "2)level=0...5"
426  << G4endl;
427 
428  fFileInitialized = true;
429 }
430 
432 {
433  return Instance()->fActiveChemistry;
434 }
435 
437 {
438  Instance()->fActiveChemistry = flag;
439 }
440 
442 {
443  return fActiveChemistry;
444 }
445 
447 {
448  fActiveChemistry = flag;
449 }
450 
452  ios_base::openmode mode)
453 {
454  if (fVerbose)
455  {
456  G4cout << "G4DNAChemistryManager: Write chemical stage into "
457  << output.data() << G4endl;
458  }
459 
460  fpgOutput_tl = new std::ofstream();
461  fpgOutput_tl->open(output.data(), mode);
462  fWriteFile = true;
463  fFileInitialized = false;
464 }
465 
467 {
468  if (fWriteFile)
469  {
470  *fpgOutput_tl << G4endl;
471  }
472 }
473 
475 {
476  if (fpgOutput_tl == 0) return;
477 
478  if (fpgOutput_tl->is_open())
479  {
480  if (fVerbose)
481  {
482  G4cout << "G4DNAChemistryManager: Close File" << G4endl;
483  }
484  fpgOutput_tl->close();
485  }
486 }
487 
490 {
491  if (!fpExcitationLevel)
492  {
494  }
495  return fpExcitationLevel;
496 }
497 
500 {
501  if (!fpIonisationLevel)
502  {
504  }
505  return fpIonisationLevel;
506 }
507 
509  G4int electronicLevel,
510  const G4Track* theIncomingTrack)
511 {
512  if (fWriteFile)
513  {
515 
516  G4double energy = -1.;
517 
518  switch (modification)
519  {
521  energy = 0;
522  break;
523  case eExcitedMolecule:
524  energy = GetExcitationLevel()->ExcitationEnergy(electronicLevel);
525  break;
526  case eIonizedMolecule:
527  energy = GetIonisationLevel()->IonisationEnergy(electronicLevel);
528  break;
529  }
530 
531  *fpgOutput_tl << setw(11) << left << theIncomingTrack->GetTrackID()
532  << setw(10) << "H2O" << left << modification << internal
533  << ":" << right << electronicLevel << left << setw(11) << ""
534  << std::setprecision(2) << std::fixed << setw(13)
535  << energy / eV << std::setprecision(6) << std::scientific
536  << setw(22)
537  << (theIncomingTrack->GetPosition().x()) / nanometer
538  << setw(22)
539  << (theIncomingTrack->GetPosition().y()) / nanometer
540  << setw(22)
541  << (theIncomingTrack->GetPosition().z()) / nanometer
542  << G4endl;
543  }
544 
545  if(fActiveChemistry)
546  {
547  G4Molecule * H2O = new G4Molecule (G4H2O::Definition());
548 
549  switch (modification)
550  {
552  H2O -> AddElectron(5,1);
553  break;
554  case eExcitedMolecule :
555  H2O -> ExciteMolecule(electronicLevel);
556  break;
557  case eIonizedMolecule :
558  H2O -> IonizeMolecule(electronicLevel);
559  break;
560  }
561 
562  G4Track * H2OTrack = H2O->BuildTrack(1*picosecond,
563  theIncomingTrack->GetPosition());
564 
565  H2OTrack -> SetParentID(theIncomingTrack->GetTrackID());
566  H2OTrack -> SetTrackStatus(fStopButAlive);
567  H2OTrack -> SetKineticEnergy(0.);
568  G4VITTrackHolder::Instance()->Push(H2OTrack);
569  }
570 // else
571 // abort();
572 }
573 
575  G4ThreeVector* finalPosition)
576 // finalPosition is a pointer because this argument is optional
577 {
578  if (fWriteFile)
579  {
581 
582  *fpgOutput_tl << setw(11) << theIncomingTrack->GetTrackID() << setw(10)
583  << "e_aq" << setw(14) << -1 << std::setprecision(2)
584  << std::fixed << setw(13)
585  << theIncomingTrack->GetKineticEnergy() / eV
586  << std::setprecision(6) << std::scientific << setw(22)
587  << (theIncomingTrack->GetPosition().x()) / nanometer
588  << setw(22)
589  << (theIncomingTrack->GetPosition().y()) / nanometer
590  << setw(22)
591  << (theIncomingTrack->GetPosition().z()) / nanometer;
592 
593  if (finalPosition != 0)
594  {
595  *fpgOutput_tl << setw(14) << (finalPosition->x()) / nanometer << setw(14)
596  << (finalPosition->y()) / nanometer << setw(14)
597  << (finalPosition->z()) / nanometer;
598  }
599 
600  *fpgOutput_tl << G4endl;
601  }
602 
603  if(fActiveChemistry)
604  {
606  G4Track * e_aqTrack(0);
607  if(finalPosition)
608  {
609  e_aqTrack = e_aq->BuildTrack(picosecond,*finalPosition);
610  }
611  else
612  {
613  e_aqTrack = e_aq->BuildTrack(picosecond,theIncomingTrack->GetPosition());
614  }
615  e_aqTrack -> SetTrackStatus(fAlive);
616  e_aqTrack -> SetParentID(theIncomingTrack->GetTrackID());
617  G4VITTrackHolder::Instance()->Push(e_aqTrack);
618  }
619 }
620 
622  double time,
623  const G4ThreeVector& position,
624  int parentID)
625 {
626  if (fWriteFile)
627  {
629 
630  *fpgOutput_tl << setw(11) << parentID << setw(10) << molecule->GetName()
631  << setw(14) << -1 << std::setprecision(2) << std::fixed
632  << setw(13) << -1 << std::setprecision(6) << std::scientific
633  << setw(22) << (position.x()) / nanometer << setw(22)
634  << (position.y()) / nanometer << setw(22)
635  << (position.z()) / nanometer;
636  *fpgOutput_tl << G4endl;
637  }
638 
639  if(fActiveChemistry)
640  {
641  G4Track* track = molecule->BuildTrack(time,position);
642  track -> SetTrackStatus(fAlive);
643  track -> SetParentID(parentID);
645  }
646  else
647  {
648  delete molecule;
649  molecule = 0;
650  }
651 }
652 
654  const G4Track* theIncomingTrack)
655 {
656  if (fWriteFile)
657  {
659 
660  *fpgOutput_tl << setw(11) << theIncomingTrack->GetTrackID() << setw(10)
661  << molecule->GetName() << setw(14) << -1
662  << std::setprecision(2) << std::fixed << setw(13)
663  << theIncomingTrack->GetKineticEnergy() / eV
664  << std::setprecision(6) << std::scientific << setw(22)
665  << (theIncomingTrack->GetPosition().x()) / nanometer
666  << setw(22)
667  << (theIncomingTrack->GetPosition().y()) / nanometer
668  << setw(22)
669  << (theIncomingTrack->GetPosition().z()) / nanometer;
670  *fpgOutput_tl << G4endl;
671  }
672 
673  if(fActiveChemistry)
674  {
675  G4Track* track = molecule->BuildTrack(theIncomingTrack->GetGlobalTime(),
676  theIncomingTrack->GetPosition());
677  track -> SetTrackStatus(fAlive);
678  track -> SetParentID(theIncomingTrack->GetTrackID());
680  }
681  else
682  {
683  delete molecule;
684  molecule = 0;
685  }
686 }
void WriteInto(const G4String &, std::ios_base::openmode mode=std::ios_base::out)
Tells the chemMan to write into a file the position and electronic state of the water molecule and th...
void Process()
Definition: G4Scheduler.cc:362
static G4VITTrackHolder * Instance()
void PushMoleculeAtParentTimeAndPlace(G4Molecule *&molecule, const G4Track *)
WARNING : In case chemistry is not activated, PushMoleculeAtParentTimeAndPlace will take care of dele...
virtual void ConstructReactionTable(G4DNAMolecularReactionTable *reactionTable)=0
static G4Electron_aq * Definition()
static G4H2O * Definition()
Definition: G4H2O.cc:46
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
CLHEP::Hep3Vector G4ThreeVector
virtual bool IsPhysicsConstructor()
void CreateSolvatedElectron(const G4Track *, G4ThreeVector *finalPosition=0)
On the same idea as the previous method but for solvated electron.
static const double nanometer
Definition: G4SIunits.hh:91
WARNING: THIS CLASS IS A PROTOTYPE G4DNAChemistryManager is called from the physics models...
const G4ThreeVector & GetPosition() const
virtual void Push(G4Track *)
G4DNAWaterIonisationStructure * fpIonisationLevel
ElectronicModification
#define G4ThreadLocal
Definition: tls.hh:84
int G4int
Definition: G4Types.hh:78
static G4Scheduler * Instance()
Definition: G4Scheduler.cc:118
static G4bool GetNewBoolValue(const char *paramString)
static G4DNAMolecularReactionTable * GetReactionTable()
G4DNAWaterExcitationStructure * GetExcitationLevel()
const G4String & GetName() const
Returns the name of the molecule.
Definition: G4Molecule.cc:303
G4Mutex chemManExistence
static G4DNAChemistryManager * fgInstance
G4double GetKineticEnergy() const
#define position
Definition: xmlparse.cc:605
G4GLOB_DLL std::ostream G4cout
virtual G4bool Notify(G4ApplicationState requestedState)
virtual void SetNewValue(G4UIcommand *, G4String)
G4VUserChemistryList * fpUserChemistryList
void PushMolecule(G4Molecule *&molecule, double time, const G4ThreeVector &position, int parentID)
WARNING : In case chemistry is not activated, PushMolecule will take care of deleting the transfered ...
bool G4bool
Definition: G4Types.hh:79
void SetGun(G4ITGun *)
Definition: G4Scheduler.hh:402
virtual void ConstructDissociationChannels()
static void DeleteInstance()
You should rather use DeleteInstance than the destructor of this class.
G4int GetTrackID() const
G4double GetGlobalTime() const
static G4GeometryManager * GetInstance()
G4bool IsMultithreadedApplication()
Definition: G4Threading.cc:134
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4bool IsWorkerThread()
Definition: G4Threading.cc:128
G4Track * BuildTrack(G4double globalTime, const G4ThreeVector &Position)
Definition: G4Molecule.cc:328
const char * data() const
static G4DNAChemistryManager * Instance()
static const double eV
Definition: G4SIunits.hh:194
G4UIcmdWithoutParameter * fpRunChem
G4int G4Mutex
Definition: G4Threading.hh:161
G4DNAWaterExcitationStructure * fpExcitationLevel
G4double energy(const ThreeVector &p, const G4double m)
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
Method used by DNA physics model to create a water molecule.
void CloseFile()
Close the file specified with WriteInto.
static void InitializeInstance()
#define G4endl
Definition: G4ios.hh:61
static void Activated(G4bool flag=true)
void OpenGeometry(G4VPhysicalVolume *vol=0)
Class Description The dynamic molecule holds all the data that change for a molecule It has a pointer...
Definition: G4Molecule.hh:93
void Gun(G4ITGun *, bool physicsTableToBuild=true)
static const double picosecond
Definition: G4SIunits.hh:142
double G4double
Definition: G4Types.hh:76
virtual void ConstructTimeStepModel(G4DNAMolecularReactionTable *reactionTable)=0
static G4DNAChemistryManager * GetInstanceIfExists()
void Initialize()
Definition: G4Scheduler.cc:299
G4bool CloseGeometry(G4bool pOptimise=true, G4bool verbose=false, G4VPhysicalVolume *vol=0)
G4DNAWaterIonisationStructure * GetIonisationLevel()
G4ApplicationState
G4UIdirectory * fpChemDNADirectory
static G4ThreadLocal std::ofstream * fpgOutput_tl
G4UIcmdWithABool * fpActivateChem
static G4ThreadLocal G4bool * fpgThreadInitialized_tl
static void DeleteInstance()