ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MolecularConfiguration.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MolecularConfiguration.cc
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 (kara (AT) cenbg . in2p3 . fr)
28 //
29 // History:
30 // -----------
31 // 10 Oct 2011 M.Karamitros created
32 //
33 // -------------------------------------------------------------------
34 
36 #include "G4MoleculeDefinition.hh"
37 #include "G4UIcommand.hh"
38 #include "G4AllocatorList.hh"
39 #include "G4AutoLock.hh"
40 #include "G4MoleculeTable.hh"
41 #include "G4Serialize.hh"
42 #include <fstream>
43 
44 using CLHEP::m2;
45 using CLHEP::s;
46 using CLHEP::kelvin;
47 
48 using namespace std;
49 
50 #if defined ( WIN32 )
51 #define __func__ __FUNCTION__
52 #endif
53 
54 /*G4ThreadLocal*/G4double G4MolecularConfiguration::fgTemperature = 298; // 310*kelvin;
55 // 25°C, used to shoot an energy
56 
57 //______________________________________________________________________________
58 // G4MolecularConfigurationManager
61 
63 
65 
67 {
68  return GetManager()->GetNumberOfCreatedSpecies();
69 }
70 
72  double,
74  molConf)
75 {
76  return molConf->fDynDiffusionCoefficient;
77 }
78 
80  const G4String& label,
81  int charge)
82 {
83  fMoleculeDefinition = moleculeDef;
84 
85  fLabel = new G4String(label);
86 
87  fMoleculeID = GetManager()->Insert(moleculeDef,
88  label,
89  this);
90  fElectronOccupancy = 0;
91 
92  fDynCharge = charge;
93 
94  fDynMass = fMoleculeDefinition->GetMass();
95 
96  fDynDiffusionCoefficient = fMoleculeDefinition->GetDiffusionCoefficient();
97  fDynVanDerVaalsRadius = fMoleculeDefinition->GetVanDerVaalsRadius();
98  fDynDecayTime = fMoleculeDefinition->GetDecayTime();
99 
100  fName = fMoleculeDefinition->GetName();
101  fName += "^";
102  fName += G4UIcommand::ConvertToString(fDynCharge);
103 
104  fFormatedName = fMoleculeDefinition->GetFormatedName();
105  fFormatedName += "^";
106  fFormatedName += "{";
107  fFormatedName += G4UIcommand::ConvertToString(fDynCharge);
108  fFormatedName += "}";
109 
111  fIsFinalized = false;
112 }
113 
115 {
116  if(fIsFinalized)
117  {
118  G4ExceptionDescription errMsg;
119  errMsg << "This molecular configuration " << GetName()
120  << " is already finalized. Therefore its "
121  " properties cannot be changed.";
122  G4Exception("G4MolecularConfiguration::MakeExceptionIfFinalized",
123  "CONF_FINALIZED",FatalException,errMsg);
124  }
125 }
126 
127 //______________________________________________________________________________
128 
131 {
132  if (!fgManager)
133  {
135  if (!fgManager) // double check for MT
136  {
137  fgManager = new G4MolecularConfiguration::
139  }
140  lock.unlock();
141  }
142 
143  return fgManager;
144 }
145 
146 //______________________________________________________________________________
147 
150 {
151 // G4cout << "Does G4AllocatorList exists= ";
152 // G4cout << (G4AllocatorList::GetAllocatorListIfExist() ? "true":"false")
153 // << G4endl;
154 
155  G4MolecularConfigurationManager::MolElectronConfTable::iterator it1;
156  G4MolecularConfigurationManager::ElectronOccupancyTable::
157  iterator it2;
158 
159  for (it1 = fElecOccTable.begin(); it1 != fElecOccTable.end(); it1++)
160  {
161  for (it2 = it1->second.begin(); it2 != it1->second.end(); it2++)
162  {
163  if (it2->second)
164  {
165  delete it2->second;
166  }
167  }
168  }
169  fElecOccTable.clear();
170  fgManager = 0;
171 }
172 
173 //______________________________________________________________________________
174 // G4MolecularConfigurationManager
178  const G4ElectronOccupancy& eOcc,
179  G4MolecularConfiguration* molConf)
180 {
181  //G4AutoLock lock(&fMoleculeCreationMutex);
182 
183  ElectronOccupancyTable& table2 = fElecOccTable[molDef];
184  ElectronOccupancyTable::iterator it = table2.find(eOcc);
185 
186  if(it == table2.end())
187  {
188  table2[eOcc] = molConf;
189  }
190  else
191  {
192  G4ExceptionDescription errMsg;
193  errMsg << "The same molecular configuration seemed to be recorded twice";
194  G4Exception("G4MolecularConfigurationManager::"
195  "SetMolecularConfiguration(const G4MoleculeDefinition* molDef,"
196  "const G4ElectronOccupancy& eOcc,"
197  "G4MolecularConfiguration* molConf)",
198  "",
200  errMsg
201  );
202  }
203 
204  fLastMoleculeID++;
205 
206  fMolConfPerID.push_back(molConf);
207 
208  //lock.unlock();
209  return fLastMoleculeID;
210 }
211 
212 //______________________________________________________________________________
213 
214 const G4ElectronOccupancy*
217  const G4ElectronOccupancy& eOcc)
218 {
219  //G4AutoLock lock(&fMoleculeCreationMutex);
220 
221  MolElectronConfTable::iterator it1 = fElecOccTable.find(molDef);
222 
223  if(it1 == fElecOccTable.end())
224  {
225  // TODO = handle exception ?
226  return 0;
227  }
228 
229  ElectronOccupancyTable& table2 = it1->second;
230  ElectronOccupancyTable::iterator it2 = table2.find(eOcc);
231 
232  //lock.unlock();
233 
234  if (it2 == table2.end())
235  {
236  // TODO = handle exception ?
237  return 0;
238  }
239 
240  return &(it2->first);
241 }
242 
243 //______________________________________________________________________________
244 
248  const G4ElectronOccupancy& eOcc)
249 {
250  MolElectronConfTable::iterator it1 = fElecOccTable.find(molDef);
251 
252  if(it1 == fElecOccTable.end()) return 0;
253 
254  ElectronOccupancyTable& table2 = it1->second;
255  ElectronOccupancyTable::iterator it = table2.find(eOcc);
256 
257  if(it == table2.end())
258  {
259  return 0;
260  }
261  else
262  {
263  return it->second;
264  }
265 
266  return 0;
267 }
268 
269 //______________________________________________________________________________
270 
273  int charge,
274  G4MolecularConfiguration* molConf)
275 {
276 
277  //G4AutoLock lock(&fMoleculeCreationMutex);
278  ChargeTable& table2 = fChargeTable[molDef];
279  ChargeTable::iterator it = table2.find(charge);
280 
281  if(it == table2.end())
282  {
283  table2[charge] = molConf;
284  }
285  else
286  {
287  //lock.unlock();
288  G4ExceptionDescription errMsg;
289  errMsg << "The same molecular configuration seemed to be recorded twice";
290  G4Exception("G4MolecularConfigurationManager::"
291  "SetMolecularConfiguration(const G4MoleculeDefinition* molDef,"
292  "int charge,"
293  "G4MolecularConfiguration* molConf)",
294  "", FatalException, errMsg);
295  }
296 
297  fLastMoleculeID++;
298  fMolConfPerID.push_back(molConf);
299  //lock.unlock();
300  return fLastMoleculeID;
301 }
302 
303 //______________________________________________________________________________
304 
308  int charge)
309 {
310  //G4AutoLock lock(&fMoleculeCreationMutex);
311 
312  MolChargeConfTable::iterator it1 = fChargeTable.find(molDef);
313 
314  if(it1 == fChargeTable.end()) return 0;
315 
316  ChargeTable& table2 = it1->second;
317  ChargeTable::iterator it = table2.find(charge);
318 
319  if(it == table2.end())
320  {
321  return 0;
322  }
323  else
324  {
325  return it->second;
326  }
327 
328  return 0;
329 
330  //lock.unlock();
331  return 0;
332 }
333 
334 //______________________________________________________________________________
335 // Static method in G4MolecularConfiguration
339 {
340  if (molDef->GetGroundStateElectronOccupancy())
341  {
342  const G4ElectronOccupancy& elecOcc =
344  G4MolecularConfiguration* molConf =
345  GetManager()->GetMolecularConfiguration(molDef, elecOcc);
346 
347  if (molConf)
348  {
349  return molConf;
350  }
351  else
352  {
353  G4MolecularConfiguration* newConf =
354  new G4MolecularConfiguration(molDef,
355  elecOcc);
356  newConf->SetUserID(molDef->GetName());
357  return newConf;
358  }
359  }
360  else
361  {
362  G4MolecularConfiguration* molConf =
363  GetManager()->GetMolecularConfiguration(molDef, molDef->GetCharge());
364  if(molConf)
365  {
366  return molConf;
367  }
368  else
369  {
370  G4MolecularConfiguration* newConf =
371  new G4MolecularConfiguration(molDef, molDef->GetCharge());
372  newConf->SetUserID(molDef->GetName());
373  return newConf;
374  }
375  }
376 }
377 
378 //______________________________________________________________________________
379 
383  const G4ElectronOccupancy& elecOcc)
384 {
385  return GetManager()->GetOrCreateMolecularConfiguration(molDef, elecOcc);
386 
387 // G4MolecularConfiguration* molConf =
388 // GetManager()->GetMolecularConfiguration(molDef, elecOcc);
389 //
390 // if (molConf)
391 // {
392 // return molConf;
393 // }
394 // else
395 // {
396 // G4MolecularConfiguration* newConf =
397 // new G4MolecularConfiguration(molDef, elecOcc);
398 // return newConf;
399 // }
400 }
401 
402 //______________________________________________________________________________
403 
407  int charge)
408 {
409  G4MolecularConfiguration* molConf =
410  GetManager()->GetMolecularConfiguration(molDef, charge);
411 
412  if(molConf)
413  {
414  return molConf;
415  }
416  else
417  {
418  G4MolecularConfiguration* newConf =
419  new G4MolecularConfiguration(molDef, charge);
420  return newConf;
421  }
422 }
423 
424 //______________________________________________________________________________
425 
427 {
429  if (fgManager) delete fgManager;
430  fgManager = 0;
431  lock.unlock();
432 }
433 
434 //______________________________________________________________________________
435 // G4MolecularConfiguration
438  const G4ElectronOccupancy& elecOcc,
439  const G4String& label)
440 {
441  fMoleculeDefinition = moleculeDef;
442 
443  fMoleculeID = GetManager()->Insert(moleculeDef,
444  elecOcc,
445  this);
446  fElectronOccupancy = GetManager()->FindCommonElectronOccupancy(moleculeDef,
447  elecOcc);
448 
449  /*
450  fgManager->fTable[fMoleculeDefinition][elecOcc] = this;
451  std::map<G4ElectronOccupancy, G4MolecularConfiguration*, comparator>::iterator it ;
452  it = fgManager->fTable[moleculeDef].find(elecOcc);
453  fElectronOccupancy = &(it->first);
454  */
455 
456  fDynCharge = fMoleculeDefinition->GetNbElectrons()
457  - fElectronOccupancy->GetTotalOccupancy()
458  + moleculeDef->GetCharge();
459  fDynMass = fMoleculeDefinition->GetMass();
460 
461  fDynDiffusionCoefficient = fMoleculeDefinition->GetDiffusionCoefficient();
462  fDynVanDerVaalsRadius = fMoleculeDefinition->GetVanDerVaalsRadius();
463  fDynDecayTime = fMoleculeDefinition->GetDecayTime();
464 
465  fName = fMoleculeDefinition->GetName();
466  fName += "^";
467  fName += G4UIcommand::ConvertToString(fDynCharge);
468 
469  fFormatedName = fMoleculeDefinition->GetFormatedName();
470  fFormatedName += "^";
471  fFormatedName += "{";
472  fFormatedName += G4UIcommand::ConvertToString(fDynCharge);
473  fFormatedName += "}";
474 
475  fLabel = 0; // let it here
476 
477  if(label != "")
478  {
479  SetLabel(label);
480  }
481 
483 
484  fIsFinalized = false;
485 }
486 
487 //______________________________________________________________________________
488 
491  int charge)
492 {
493  fMoleculeDefinition = moleculeDef;
494 
495  fMoleculeID = GetManager()->Insert(moleculeDef,
496  charge,
497  this);
498  fElectronOccupancy = 0;
499 
500  fDynCharge = charge;
501  fDynMass = fMoleculeDefinition->GetMass();
502 
503  fDynDiffusionCoefficient = fMoleculeDefinition->GetDiffusionCoefficient();
504  fDynVanDerVaalsRadius = fMoleculeDefinition->GetVanDerVaalsRadius();
505  fDynDecayTime = fMoleculeDefinition->GetDecayTime();
506 
507  fName = fMoleculeDefinition->GetName();
508  fName += "^";
509  fName += G4UIcommand::ConvertToString(fDynCharge);
510 
511  fFormatedName = fMoleculeDefinition->GetFormatedName();
512  fFormatedName += "^";
513  fFormatedName += "{";
514  fFormatedName += G4UIcommand::ConvertToString(fDynCharge);
515  fFormatedName += "}";
516 
517  fLabel = 0;
518 
520 
521  fIsFinalized = false;
522 }
523 
524 //______________________________________________________________________________
525 
527 {
528  if (fgManager) fgManager->RemoveMolecularConfigurationFromTable(this);
529 
530 // if (G4AllocatorList::GetAllocatorListIfExist())
531 // {
532 // if (fElectronOccupancy)
533 // {
534 // delete fElectronOccupancy;
535 // fElectronOccupancy = 0;
536 // }
537 // }
538 }
539 
540 //______________________________________________________________________________
541 
544 ChangeConfiguration(const G4ElectronOccupancy& newElectronOccupancy) const
545 {
546  G4MolecularConfiguration* output =
547  GetManager()->GetMolecularConfiguration(fMoleculeDefinition,
548  newElectronOccupancy);
549 
550  if (!output)
551  {
552  output = new G4MolecularConfiguration(fMoleculeDefinition,
553  newElectronOccupancy);
554  }
555  return output;
556 }
557 
558 //______________________________________________________________________________
559 
562 {
563  G4MolecularConfiguration* output =
564  GetManager()->GetMolecularConfiguration(fMoleculeDefinition, charge);
565 
566  if (!output)
567  {
568  output = new G4MolecularConfiguration(fMoleculeDefinition, charge);
569  }
570  return output;
571 }
572 
573 //______________________________________________________________________________
574 
577 {
578 // if (&right == this) return *this;
579  return *this;
580 }
581 
582 //______________________________________________________________________________
583 
588 {
589 // MakeExceptionIfFinalized();
590  CheckElectronOccupancy(__func__);
591  G4ElectronOccupancy newElectronOccupancy(*fElectronOccupancy);
592 
593  newElectronOccupancy.RemoveElectron(ExcitedLevel, 1);
594  newElectronOccupancy.AddElectron(5, 1);
595 
596  return ChangeConfiguration(newElectronOccupancy);
597 }
598 
599 //______________________________________________________________________________
600 
605 {
606 // MakeExceptionIfFinalized();
607  CheckElectronOccupancy(__func__);
608  G4ElectronOccupancy newElectronOccupancy(*fElectronOccupancy);
609 
610  if (newElectronOccupancy.GetOccupancy(IonizedLevel) != 0)
611  {
612  newElectronOccupancy.RemoveElectron(IonizedLevel, 1);
613  }
614  else
615  {
616  G4String errMsg = "There is no electron on the orbit "
617  + G4UIcommand::ConvertToString(IonizedLevel)
618  + " you want to free. The molecule's name you want to ionized is "
619  + GetName();
620  G4Exception("G4MolecularConfiguration::IonizeMolecule",
621  "",
623  errMsg);
624  PrintState();
625  }
626 
627  // DEBUG
628  // PrintState();
629 
630  return ChangeConfiguration(newElectronOccupancy);
631 }
632 
633 //______________________________________________________________________________
634 
636  G4int number) const
637 {
638 // MakeExceptionIfFinalized();
639  CheckElectronOccupancy(__func__);
640  G4ElectronOccupancy newElectronOccupancy(*fElectronOccupancy);
641  newElectronOccupancy.AddElectron(orbit, number);
642  return ChangeConfiguration(newElectronOccupancy);
643 }
644 
645 //______________________________________________________________________________
646 
649  G4int number) const
650 {
651 // MakeExceptionIfFinalized();
652  CheckElectronOccupancy(__func__);
653  G4ElectronOccupancy newElectronOccupancy(*fElectronOccupancy);
654 
655  if (newElectronOccupancy.GetOccupancy(orbit) != 0)
656  {
657  newElectronOccupancy.RemoveElectron(orbit, number);
658  }
659  else
660  {
661  G4String errMsg = "There is already no electron into the orbit "
663  + " you want to free. The molecule's name is " + GetName();
664  G4Exception("G4MolecularConfiguration::RemoveElectron",
665  "",
666  JustWarning,
667  errMsg);
668  PrintState();
669  }
670 
671  return ChangeConfiguration(newElectronOccupancy);
672 }
673 
674 //______________________________________________________________________________
675 
678  G4int orbitToFill) const
679 {
680 // MakeExceptionIfFinalized();
681  CheckElectronOccupancy(__func__);
682  G4ElectronOccupancy newElectronOccupancy(*fElectronOccupancy);
683 
684  if (newElectronOccupancy.GetOccupancy(orbitToFree) >= 1)
685  {
686  newElectronOccupancy.RemoveElectron(orbitToFree, 1);
687  newElectronOccupancy.AddElectron(orbitToFill, 1);
688  }
689  else
690  {
691  G4String errMsg = "There is no electron on the orbit "
692  + G4UIcommand::ConvertToString(orbitToFree)
693  + " you want to free. The molecule's name is " + GetName();
694  G4Exception("G4MolecularConfiguration::MoveOneElectron",
695  "",
697  errMsg);
698  PrintState();
699  }
700 
701  return ChangeConfiguration(newElectronOccupancy);
702 }
703 
704 //______________________________________________________________________________
705 
707 {
708 // if (fName.isNull())
709 // {
710 // fName = fMoleculeDefinition->GetName();
711 // fName += "^";
712 // // fName+= "{";
713 // fName += G4UIcommand::ConvertToString(fDynCharge);
714 // // fName+= "}";
715 // }
716  return fName;
717 }
718 
719 //______________________________________________________________________________
720 
722 {
723 // if (fFormatedName.isNull())
724 // {
725 // fFormatedName = fMoleculeDefinition->GetFormatedName();
726 // fFormatedName += "^";
727 // fFormatedName += "{";
728 // fFormatedName += G4UIcommand::ConvertToString(fDynCharge);
729 // fFormatedName += "}";
730 // }
731  return fFormatedName;
732 }
733 
734 //______________________________________________________________________________
735 
737 {
738  return fMoleculeDefinition->GetAtomsNumber();
739 }
740 
741 //______________________________________________________________________________
742 
744 {
745  CheckElectronOccupancy(__func__);
746  return fElectronOccupancy->GetTotalOccupancy();
747 }
748 
749 //______________________________________________________________________________
750 
752 {
753  G4cout << "-------------- Start Printing State " << GetName()
754  << " ---------------" << G4endl;
755 
756  if (fElectronOccupancy)
757  {
758  G4cout << "--------------Print electronic state of " << GetName()
759  << "---------------" << G4endl;
760  fElectronOccupancy->DumpInfo();
761  if(fElectronOccupancy==fMoleculeDefinition->GetGroundStateElectronOccupancy())
762  {
763  G4cout<<"At ground state"<<G4endl;
764  }
765  }
766  else
767  {
768  G4cout << "--- No electron occupancy set up ---" << G4endl;
769  }
770 
771  G4cout << "Charge :"
772  << fDynCharge
773  << G4endl;
774 
775  if(fLabel)
776  {
777  G4cout << "Label :"
778  << GetLabel()
779  << G4endl;
780  }
781  G4cout << "-------------- End Of State " << GetName()
782  << " -----------------------" << G4endl;
783 }
784 
785 //______________________________________________________________________________
786 
787 // added - to be transformed in a "Decay method"
788 const vector<const G4MolecularDissociationChannel*>*
790 {
791  // if (fElectronOccupancy == 0) return 0;
792  return fMoleculeDefinition->GetDecayChannels(this);
793 }
794 
795 //______________________________________________________________________________
796 
798 {
799  if(fMoleculeDefinition) return fMoleculeDefinition->GetPDGEncoding();
800  else G4Exception("G4MolecularConfiguration::GetMoleculeID",
801  "",
803  "You should first enter a molecule definition");
804 
805  return INT_MAX;
806 }
807 
808 //______________________________________________________________________________
809 
810 const char* removePath(const char* path)
811 {
812  const char* pDelimeter = strrchr(path, '\\');
813  if (pDelimeter) path = pDelimeter + 1;
814 
815  pDelimeter = strrchr(path, '/');
816  if (pDelimeter) path = pDelimeter + 1;
817 
818  return path;
819 }
820 
821 //______________________________________________________________________________
822 
823 void G4MolecularConfiguration::CheckElectronOccupancy(const char* function) const
824 {
825  if (fElectronOccupancy == 0)
826  {
827  G4String functionName(function);
828  G4ExceptionDescription description;
829  description
830  << "No G4ElectronOccupancy was defined for molecule definition : "
831  << fMoleculeDefinition->GetName()
832  << ". The definition was probably defined using the charge state, "
833  "rather than electron state.";
834 
835  G4Exception(functionName, "", FatalErrorInArgument, description);
836  }
837 }
838 
839 //______________________________________________________________________________
840 
843 {
844  //G4AutoLock lock(&fMoleculeCreationMutex);
845 
846  LabelTable& tmpMap = fLabelTable[molConf->fMoleculeDefinition];
847 
848  LabelTable::iterator it = tmpMap.find(*molConf->fLabel);
849 
850  if(it == tmpMap.end())
851  {
852  tmpMap[*(molConf->fLabel)] = molConf;
853  }
854  else
855  {
856  G4ExceptionDescription errMsg;
857  errMsg << "The same molecular configuration seemed to be recorded twice";
858  G4Exception("G4MolecularConfigurationManager::"
859  "SetMolecularConfiguration(const G4MoleculeDefinition* molDef,"
860  "const G4String& label,"
861  "G4MolecularConfiguration* molConf)",
862  "", FatalException, errMsg);
863  }
864 
865  //lock.unlock();
866 }
867 
870 {
871  UserIDTable::iterator it = fUserIDTable.find(userID);
872 
873  if(it == fUserIDTable.end())
874  {
875  fUserIDTable[userID] = molecule;
876  }
877  else if(molecule != it->second)
878  {
879  // TODO improve exception
880  // exception
881  G4ExceptionDescription description;
882  description << "The user identifier " << userID
883  << " was already given in another configuration in the table"
884  << G4endl;
885  G4Exception("G4MolecularConfiguration::G4MolecularConfigurationManager::AddUserID",
886  "CONF_ALREADY_RECORDED",
888  description);
889  }
890 }
891 
892 //______________________________________________________________________________
893 
896 {
897  MolElectronConfTable::iterator it1 =
898  fElecOccTable.find(configuration->GetDefinition());
899  MolElectronConfTable::iterator end = fElecOccTable.end();
900 
901  if (it1 == end) return;
902 
903  std::map<G4ElectronOccupancy, G4MolecularConfiguration*, comparator>::
904  iterator it2 =
905  it1->second.find(*configuration->GetElectronOccupancy());
906 
907  if (it2 == it1->second.end()) return;
908 
909  it2->second = 0;
910 // it1->second.erase(it2);
911 
912  configuration->fElectronOccupancy = 0;
913 }
914 
915 //______________________________________________________________________________
916 
920  const G4String& label)
921 {
922  //G4AutoLock lock(&fMoleculeCreationMutex);
923 
924  MolLabelConfTable::iterator it1 = fLabelTable.find(molDef);
925 
926  if(it1 == fLabelTable.end()) return 0;
927 
928  LabelTable& table2 = it1->second;
929 
930  LabelTable::iterator it2 = table2.find(label);
931 
932  //lock.unlock();
933 
934  if(it2 == table2.end()) return 0;
935  return it2->second;
936 }
937 
938 //______________________________________________________________________________
939 
943 {
944  if(moleculeID > (int) fMolConfPerID.size() ||
945  moleculeID < 0) return 0;
946 
947  return fMolConfPerID[moleculeID];
948 }
949 
950 //______________________________________________________________________________
951 
952 G4int
955  const G4String& label,
956  G4MolecularConfiguration* molConf)
957 {
958  G4AutoLock lock(&fMoleculeCreationMutex);
959  LabelTable& tmpMap = fLabelTable[molDef];
960  LabelTable::iterator it = tmpMap.find(label);
961 
962  if(it == tmpMap.end())
963  {
964  fLastMoleculeID++;
965  tmpMap[label] = molConf;
966  lock.unlock();
967  }
968  else
969  {
970  lock.unlock();
971  G4ExceptionDescription errMsg;
972  errMsg << "The same molecular configuration seemed to be recorded twice";
973  G4Exception("G4MolecularConfigurationManager::"
974  "SetMolecularConfiguration(const G4MoleculeDefinition* molDef,"
975  "const G4String& label,"
976  "G4MolecularConfiguration* molConf)",
977  "", FatalException, errMsg);
978  }
979 
980  fMolConfPerID.push_back(molConf);
981 
982  return fLastMoleculeID;
983 }
984 
985 //______________________________________________________________________________
986 
989  const G4String& label)
990 {
991  return GetManager()->GetMolecularConfiguration(molDef, label);
992 }
993 
994 //______________________________________________________________________________
995 
998 {
999  return GetManager()->GetMolecularConfiguration(moleculeID);
1000 }
1001 
1002 //______________________________________________________________________________
1003 
1006  const G4MoleculeDefinition* molDef,
1007  int charge,
1008  const G4String& label,
1009  bool& wasAlreadyCreated)
1010 {
1011  wasAlreadyCreated = false;
1012  G4MolecularConfiguration* molConf =
1013  GetManager()->GetMolecularConfiguration(molDef, charge);
1014 
1015  if (molConf)
1016  {
1017  if(molConf->fLabel == 0)
1018  {
1019  molConf->SetLabel(label);
1020  G4ExceptionDescription wMsg ;
1021  wMsg << "The molecular configuration for the definition named "
1022  << molDef->GetName()
1023  << " with charge " << charge << " has already been created "
1024  "but with NO label";
1025  G4Exception("G4MolecularConfiguration::CreateMolecularConfiguration",
1026  "DOUBLE_CREATION",
1027  JustWarning,
1028  wMsg);
1029  }
1030  else if(*(molConf->fLabel) == "" )
1031  {
1032  molConf->SetLabel(label);
1033  }
1034  else if(*(molConf->fLabel) != label)
1035  {
1036  G4ExceptionDescription errMsg ;
1037  errMsg << "The molecular configuration for the definition named "
1038  << molDef->GetName()
1039  << " with charge " << charge << " has already been created "
1040  "but with a different label :"
1041  << molConf->GetLabel();
1042  G4Exception("G4MolecularConfiguration::CreateMolecularConfiguration",
1043  "DOUBLE_CREATION",
1045  errMsg);
1046  // KILL APP
1047  }
1048 
1049  if(molConf->fUserIdentifier == "")
1050  {
1051  molConf->fUserIdentifier = userIdentifier;
1052 
1053  G4ExceptionDescription wMsg ;
1054  wMsg << "The molecular configuration for the definition named "
1055  << molDef->GetName()
1056  << " with label " << label << " has already been created.";
1057  G4Exception("G4MolecularConfiguration::CreateMolecularConfiguration",
1058  "DOUBLE_CREATION",
1059  JustWarning,
1060  wMsg);
1061  }
1062  else if(molConf->fUserIdentifier != userIdentifier)
1063  {
1064  G4ExceptionDescription errMsg ;
1065  errMsg << "The molecular configuration for the definition named "
1066  << molDef->GetName()
1067  << " with label " << label << " has already been created "
1068  "BUT with a different user ID :"
1069  << molConf->fUserIdentifier;
1070  G4Exception("G4MolecularConfiguration::CreateMolecularConfiguration",
1071  "DOUBLE_CREATION",
1073  errMsg);
1074  // KILL APP
1075  }
1076 
1077  wasAlreadyCreated = true;
1078  return molConf;
1079  }
1080  else
1081  {
1082  G4MolecularConfiguration* newConf =
1083  new G4MolecularConfiguration(molDef, label, charge);
1084  newConf->fUserIdentifier = userIdentifier;
1085 
1086  GetManager()->AddUserID(userIdentifier, newConf);
1087 
1088 // G4MoleculeTable::Instance()->RecordMolecularConfiguration(userIdentifier,
1089 // newConf);
1090  return newConf;
1091  }
1092 }
1093 
1094 //______________________________________________________________________________
1095 
1099  const G4MoleculeDefinition* molDef,
1100  bool& wasAlreadyCreated)
1101 {
1102  wasAlreadyCreated = false;
1103  G4MolecularConfiguration* preRegisteredMolConf =
1104  GetManager()->GetMolecularConfiguration(userIdentifier);
1105 
1106  if(preRegisteredMolConf)
1107  {
1108  if(preRegisteredMolConf->GetDefinition() == molDef)
1109  {
1110  wasAlreadyCreated = true;
1111  return preRegisteredMolConf;
1112  }
1113  }
1114 
1115  if(molDef->GetGroundStateElectronOccupancy())
1116  {
1117  const G4ElectronOccupancy& elecOcc = *molDef
1119  G4MolecularConfiguration* molConf =
1120  GetManager()->GetMolecularConfiguration(molDef, elecOcc);
1121 
1122  if(molConf)
1123  {
1124  if(molConf->fUserIdentifier == "")
1125  {
1126  molConf->fUserIdentifier = userIdentifier;
1127  }
1128  else if(molConf->fUserIdentifier != userIdentifier)
1129  {
1130  G4ExceptionDescription errMsg;
1131  errMsg << "A molecular configuration for the definition named "
1132  << molDef->GetName() << " has already been created "
1133  "and recorded with a different user ID "
1134  << molConf->fUserIdentifier;
1135  G4Exception("G4MolecularConfiguration::CreateMolecularConfiguration",
1136  "DOUBLE_CREATION",
1138  errMsg);
1139  }
1140 // TODO exception
1141  G4ExceptionDescription errMsg;
1142  errMsg << "A molecular configuration for the definition named "
1143  << molDef->GetName() << " has already been created.";
1144  G4Exception("G4MolecularConfiguration::CreateMolecularConfiguration",
1145  "DOUBLE_CREATION",
1146  JustWarning,
1147  errMsg);
1148  wasAlreadyCreated = true;
1149  return molConf;
1150  }
1151  else
1152  {
1153  // G4cout << "Create molConf for " << molDef->GetName() << G4endl;
1154  G4MolecularConfiguration* newConf = new G4MolecularConfiguration(molDef,
1155  elecOcc);
1156  newConf->fUserIdentifier = userIdentifier;
1157 
1158  GetManager()->AddUserID(userIdentifier, newConf);
1159 
1160 // G4MoleculeTable::Instance()->RecordMolecularConfiguration(userIdentifier,
1161 // newConf);
1162  return newConf;
1163  }
1164  }
1165  else
1166  {
1167  return CreateMolecularConfiguration(userIdentifier,
1168  molDef,
1169  molDef->GetName(),
1170  molDef->GetCharge(),
1171  wasAlreadyCreated);
1172  }
1173 }
1174 
1175 //______________________________________________________________________________
1176 
1180  const G4MoleculeDefinition* molDef,
1181  const G4String& label,
1182  bool& wasAlreadyCreated)
1183 {
1184  assert(label != "");
1185  wasAlreadyCreated = false;
1186 
1187  G4MolecularConfiguration* molConf =
1188  GetManager()->GetMolecularConfiguration(molDef, label);
1189  if(molConf)
1190  {
1191  if(molConf->fLabel
1192  && *molConf->fLabel == label)
1193  {
1194  wasAlreadyCreated = true;
1195  return molConf;
1196  }
1197  else if(molConf->fLabel == 0)
1198  {
1199  wasAlreadyCreated = true;
1200  molConf->SetLabel(label);
1201  return molConf;
1202  }
1203  else if(*molConf->fLabel == "")
1204  {
1205  wasAlreadyCreated = true;
1206  molConf->SetLabel(label);
1207  return molConf;
1208  }
1209 
1210  molConf->PrintState();
1211  G4ExceptionDescription errMsg ;
1212  errMsg << "A molecular configuration for the definition named "
1213  << molDef->GetName()
1214  << " has already been created "
1215  "with user ID "
1216  << molConf->fUserIdentifier << " and label "
1217  << molConf->GetLabel();
1218  G4Exception("G4MolecularConfiguration::CreateMolecularConfiguration",
1219  "DOUBLE_CREATION",
1221  errMsg);
1222  // KILL APP
1223  }
1224  else
1225  {
1226  G4MolecularConfiguration* newConf =
1227  new G4MolecularConfiguration(molDef,
1228  label,
1229  molDef->GetCharge());
1230  newConf->fUserIdentifier = userIdentifier;
1231 
1232  GetManager()->AddUserID(userIdentifier, newConf);
1233 
1234 // G4MoleculeTable::Instance()->
1235 // RecordMolecularConfiguration(userIdentifier, newConf);
1236  return newConf;
1237  }
1238  return molConf;
1239 }
1240 
1241 //______________________________________________________________________________
1242 
1246  const G4MoleculeDefinition* molDef,
1247  const G4String& label,
1248  const G4ElectronOccupancy& eOcc,
1249  bool& wasAlreadyCreated)
1250 {
1251  assert(label != "");
1252  wasAlreadyCreated = false;
1253 
1254  G4MolecularConfiguration* molConf =
1255  GetManager()->GetMolecularConfiguration(molDef, eOcc);
1256 
1257  if(molConf)
1258  {
1259  if(molConf->GetElectronOccupancy())
1260  {
1261  if(*molConf->GetElectronOccupancy() == eOcc)
1262  {
1263  if(molConf->fLabel && *molConf->fLabel == label)
1264  {
1265  wasAlreadyCreated = true;
1266  return molConf;
1267  }
1268  else if(molConf->fLabel == 0)
1269  {
1270  wasAlreadyCreated = true;
1271  molConf->SetLabel(label);
1272  return molConf;
1273  }
1274  else if(*molConf->fLabel == "")
1275  {
1276  wasAlreadyCreated = true;
1277  molConf->SetLabel(label);
1278  return molConf;
1279  }
1280  }
1281  }
1282 
1283 
1284  molConf->PrintState();
1285  G4ExceptionDescription errMsg ;
1286  errMsg << "A molecular configuration for the definition named "
1287  << molDef->GetName()
1288  << " has already been created "
1289  "with user ID "
1290  << molConf->fUserIdentifier
1291  << " and possible different electronic state";
1292  G4Exception("G4MolecularConfiguration::CreateMolecularConfiguration",
1293  "DOUBLE_CREATION",
1295  errMsg);
1296  }
1297  else
1298  {
1299  G4MolecularConfiguration* newConf =
1300  new G4MolecularConfiguration(molDef,
1301  eOcc,
1302  label);
1303  newConf->fUserIdentifier = userIdentifier;
1304 
1305  GetManager()->AddUserID(userIdentifier, newConf);
1306 
1307 // G4MoleculeTable::Instance()->
1308 // RecordMolecularConfiguration(userIdentifier, newConf);
1309  return newConf;
1310  }
1311  return molConf;
1312 }
1313 
1314 
1315 //______________________________________________________________________________
1316 
1320  const G4ElectronOccupancy& eOcc)
1321 {
1322  MolElectronConfTable::iterator it1 = fElecOccTable.find(molDef);
1323 
1324  if(it1 == fElecOccTable.end())
1325  {
1326  return new G4MolecularConfiguration(molDef, eOcc);
1327  }
1328 
1329  ElectronOccupancyTable& table2 = it1->second;
1330  ElectronOccupancyTable::iterator it = table2.find(eOcc);
1331 
1332  if(it == table2.end())
1333  {
1334  G4MolecularConfiguration* molConf =
1335  new G4MolecularConfiguration(molDef, eOcc);
1336 // molConf->Finalize();
1337  return molConf;
1338  }
1339  else
1340  {
1341  return it->second;
1342  }
1343 
1344  return 0;
1345 }
1346 
1347 //______________________________________________________________________________
1348 
1352  int charge)
1353 {
1354  MolChargeConfTable::iterator it1 = fChargeTable.find(molDef);
1355 
1356  if(it1 == fChargeTable.end())
1357  {
1358  G4AutoLock lock(&fMoleculeCreationMutex);
1359 
1360  G4MolecularConfiguration* newConf = new G4MolecularConfiguration(molDef, charge);
1361  return newConf ;
1362  }
1363 
1364  ChargeTable& table2 = it1->second;
1365  ChargeTable::iterator it = table2.find(charge);
1366 
1367  if(it == table2.end())
1368  {
1369  G4AutoLock lock(&fMoleculeCreationMutex);
1370 
1371  G4MolecularConfiguration* newConf =
1372  new G4MolecularConfiguration(molDef, charge);
1373 // newConf->Finalize();
1374  return newConf ;
1375  }
1376  else
1377  {
1378  return it->second;
1379  }
1380 
1381  return 0;
1382 }
1383 
1384 //______________________________________________________________________________
1385 
1387 {
1388  G4String moleculeName = fMoleculeDefinition->GetName();
1389  WRITE(out, moleculeName);
1390 
1391 // if(fLabel)
1392 // out << fLabel;
1393 // else
1394 // out << "";
1395  WRITE(out,fDynDiffusionCoefficient);
1396  WRITE(out,fDynVanDerVaalsRadius);
1397  WRITE(out,fDynDecayTime);
1398  WRITE(out,fDynMass);
1399  WRITE(out,fDynCharge);
1400  WRITE(out,fMoleculeID);
1401  WRITE(out,fFormatedName);
1402  WRITE(out,fName);
1403  WRITE(out,fIsFinalized);
1404 }
1405 
1406 //______________________________________________________________________________
1407 
1409 {
1410  G4String moleculeName;
1411  READ(in, moleculeName);
1412  fMoleculeDefinition =
1414 
1415 // G4String label;
1416 //
1417 // in.read((char*)(&label), sizeof(label));
1418 //
1419 // if(label)
1420 // fLabel = new G4String(label);
1421 // else
1422 // fLabel = 0;
1423  READ(in,fDynDiffusionCoefficient);
1424  READ(in,fDynVanDerVaalsRadius);
1425  READ(in,fDynDecayTime);
1426  READ(in,fDynMass);
1427  READ(in,fDynCharge);
1428  READ(in,fMoleculeID);
1429  READ(in,fFormatedName);
1430  READ(in,fName);
1431  READ(in,fIsFinalized);
1432 }
1433 
1434 //______________________________________________________________________________
1435 
1437 {
1438  return new G4MolecularConfiguration(in);
1439 }
1440 
1441 //______________________________________________________________________________
1442 
1444 {
1445  fLabel = 0; // TODO: for now not serialized
1446  Unserialize(in);
1447  fMoleculeDefinition = 0;
1448  fElectronOccupancy = 0;
1449  if(fElectronOccupancy)
1450  {
1451  GetManager()->Insert(fMoleculeDefinition, *fElectronOccupancy, this);
1452  fElectronOccupancy =
1453  GetManager()->FindCommonElectronOccupancy(fMoleculeDefinition,
1454  *fElectronOccupancy);
1455 
1456  if(fLabel)
1457  {
1458  GetManager()->RecordNewlyLabeledConfiguration(this);
1459  }
1460  }
1461  else if(fLabel)
1462  {
1463  fMoleculeID = GetManager()->Insert(fMoleculeDefinition, *fLabel, this);
1464  }
1465  else if(fDynCharge)
1466  {
1467  fMoleculeID = GetManager()->Insert(fMoleculeDefinition, fDynCharge, this);
1468  }
1469 }
1470 
1471 //______________________________________________________________________________
1472 
1474 {
1475  fUserIdentifier = userID;
1476  GetManager()->AddUserID(userID, this);
1477 // G4MoleculeTable::Instance()->RecordMolecularConfiguration(userID, this);
1478 }
1479 
1480 //______________________________________________________________________________
1481 
1482 double G4MolecularConfiguration::DiffCoeffWater(double temperature_K)
1483 {
1484  return pow(10, 4.311
1485  - 2.722e3/temperature_K
1486  + 8.565e5/(temperature_K *temperature_K)
1487  - 1.181e8/(temperature_K*temperature_K*temperature_K ))*1e-9*m2/s;
1488 }
1489 
1490 //______________________________________________________________________________
1491 
1492 void
1495 {
1496  double D_water_0 = DiffCoeffWater(fgTemperature);
1497  double D_water_f = DiffCoeffWater(temperature_K);
1498 
1499  G4cout << "Scaling factor = " << D_water_f/D_water_0 << G4endl;
1500 
1503 
1504  while(it())
1505  {
1506  G4MolecularConfiguration* conf = it.value();
1507  double D_0 = conf->GetDiffusionCoefficient() ;
1508  double D_f = D_water_f * D_0 /D_water_0;
1509  conf->SetDiffusionCoefficient(D_f);
1510  };
1511 }
1512 
1513 //______________________________________________________________________________
1514 
1516 {
1517  if(bool(fDiffParam) == false)
1518  {
1520  }
1521 }
1522 
1523 //______________________________________________________________________________
1524 
1526 {
1527  ScaleAllDiffusionCoefficientsOnWater(temperature);
1528  fgTemperature = temperature;
1529 }
1530 
1531 //______________________________________________________________________________
1532 
1534 {
1535  return fgTemperature;
1536 }
1537 
1538 //______________________________________________________________________________
1539 
1543 {
1544  for(auto it : fMolConfPerID)
1545  {
1546  if(it->GetUserID() == userID) return it;
1547  }
1548  return 0;
1549 }
1550 
1551 //______________________________________________________________________________
1552 
1555 {
1556  return GetManager()->GetMolecularConfiguration(userID);
1557 }
1558 
1559 //______________________________________________________________________________
1560 
1562 {
1563  const std::vector<G4MolecularConfiguration*>& species =
1564  GetManager()->GetAllSpecies();
1565 
1566  for(size_t i = 0; i < species.size() ; ++i)
1567  {
1568  species[i]->Finalize();
1569  }
1570 
1571 }