ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAMolecularReactionTable.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNAMolecularReactionTable.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 // WARNING : This class is released as a prototype.
30 // It might strongly evolve or even disapear in the next releases.
31 //
32 // History:
33 // -----------
34 // 10 Oct 2011 M.Karamitros created
35 //
36 // -------------------------------------------------------------------
37 
38 #include <iomanip>
39 
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4UIcommand.hh"
44 #include "G4VDNAReactionModel.hh"
46 #include "G4MoleculeTable.hh"
49 #include "G4IosFlagsSaver.hh"
50 #include "G4Exp.hh"
51 
52 using namespace std;
53 
55 
57  : fpReactant1(nullptr)
58  , fpReactant2(nullptr)
59  , fObservedReactionRate(0.)
60  , fEffectiveReactionRadius(0.)
61  , fReactionID(0)
62 {
63 }
64 
66  Reactant* pReactant1,
67  Reactant* pReactant2)
68  : fpReactant1(pReactant1)
69  , fpReactant2(pReactant2)
70  , fObservedReactionRate(reactionRate)
71  , fEffectiveReactionRadius(0.)
72  , fReactionID(0)
73 {
75 }
76 
78  const G4String& reactant1,
79  const G4String& reactant2)
80  : fpReactant1(nullptr)
81  , fpReactant2(nullptr)
82  , fObservedReactionRate(reactionRate)
83  , fEffectiveReactionRadius(0.)
84  , fReactionID(0)
85 {
86  SetReactant1(reactant1);
87  SetReactant2(reactant2);
89 }
90 
92 {
93  fProducts.clear();
94 }
95 
97 {
98  G4double sumDiffCoeff = 0.;
99 
100  if (fpReactant1 == fpReactant2)
101  {
102  sumDiffCoeff = fpReactant1->GetDiffusionCoefficient();
104  }
105  else
106  {
107  sumDiffCoeff = fpReactant1->GetDiffusionCoefficient()
110  }
111 }
112 
114 {
115  return fReactionID;
116 }
117 
119 {
120  fReactionID = ID;
121 }
122 
124 {
125  fpReactant1 = pReactive;
126 }
127 
129 {
130  fpReactant2 = pReactive;
131 }
132 
134  Reactant* pReactant2)
135 {
136  fpReactant1 = pReactant1;
137  fpReactant2 = pReactant2;
138 }
139 
141 {
142  fProducts.push_back(pMolecule);
143 }
144 
146 {
147  return fProducts.size();
148 }
149 
151 {
152  return fProducts[i];
153 }
154 
156 {
157  return &fProducts;
158 }
159 
161 {
162  fProducts.clear();
163 }
164 
166 {
168 }
170 {
172 }
174  const G4String& reactant2)
175 {
178 }
179 
181 {
182  return std::make_pair(fpReactant1, fpReactant2);
183 }
184 
186 {
187  return fpReactant1;
188 }
189 
191 {
192  return fpReactant2;
193 }
194 
196 {
197  fObservedReactionRate = rate;
198 }
199 
201 {
202  return fObservedReactionRate;
203 }
204 
206 {
208 }
209 
211 {
213 }
214 
216 {
217  fProducts.push_back(G4MoleculeTable::Instance()->GetConfiguration(molecule));
218 }
219 
220 double G4DNAMolecularReactionData::PolynomialParam(double temp_K, std::vector<double> P)
221 {
222  double inv_temp = 1. / temp_K;
223 
224  return pow(10,
225  P[0] + P[1] * inv_temp + P[2] * pow(inv_temp, 2)
226  + P[3] * pow(inv_temp, 3) + P[4] * pow(inv_temp, 4))
227  * (1e-3 * CLHEP::m3 / (CLHEP::mole * CLHEP::s));
228 }
229 
230 double G4DNAMolecularReactionData::ArrehniusParam(double temp_K, std::vector<double> P)
231 {
232  return P[0] * G4Exp(P[1] / temp_K)*
233  (1e-3 * CLHEP::m3 / (CLHEP::mole * CLHEP::s));
234 }
235 
237  double temp_init,
238  double rateCste_init)
239 {
240  double D0 = G4MolecularConfiguration::DiffCoeffWater(temp_init);
241  double Df = G4MolecularConfiguration::DiffCoeffWater(temp_K);
242  return Df * rateCste_init / D0;
243 }
244 
245 //==============================================================================
246 // REACTION TABLE
247 //==============================================================================
248 
250 {
251  if (!fpInstance)
252  {
254  }
255  return fpInstance;
256 }
257 
258 //_____________________________________________________________________________________
259 
261 {
262  if (!fpInstance)
263  {
265  }
266  return fpInstance;
267 }
268 
269 //_____________________________________________________________________________________
270 
272 {
273  if (fpInstance)
274  {
275  delete fpInstance;
276  }
277  fpInstance = nullptr;
278 }
279 
280 //_____________________________________________________________________________________
281 
284  , fVerbose(false)
285  , fpMessenger(new G4ReactionTableMessenger(this))
286 {
287 }
288 
290 
292 {
293  const auto pReactant1 = pReactionData->GetReactant1();
294  const auto pReactant2 = pReactionData->GetReactant2();
295 
296  fReactionData[pReactant1][pReactant2] = pReactionData;
297  fReactantsMV[pReactant1].push_back(pReactant2);
298  fReactionDataMV[pReactant1].push_back(pReactionData);
299 
300  if (pReactant1 != pReactant2)
301  {
302  fReactionData[pReactant2][pReactant1] = pReactionData;
303  fReactantsMV[pReactant2].push_back(pReactant1);
304  fReactionDataMV[pReactant2].push_back(pReactionData);
305  }
306 
307  fVectorOfReactionData.emplace_back(pReactionData);
308  pReactionData->SetReactionID(fVectorOfReactionData.size());
309 }
310 
311 //_____________________________________________________________________________________
312 
314  Reactant* pReactant1,
315  Reactant* pReactant2)
316 {
317  auto reactionData = new G4DNAMolecularReactionData(reactionRate, pReactant1, pReactant2);
318  SetReaction(reactionData);
319 }
320 
321 //_____________________________________________________________________________________
322 
324 {
325  // Print Reactions and Interaction radius for jump step = 3ps
326 
327  G4IosFlagsSaver iosfs(G4cout);
328 
329  if (pReactionModel && !(pReactionModel->GetReactionTable()))
330  {
331  pReactionModel->SetReactionTable(this);
332  }
333 
334  ReactivesMV::iterator itReactives;
335 
336  std::map<Reactant*, std::map<Reactant*, G4bool>> alreadyPrint;
337 
338  G4cout << "Number of chemical species involved in reactions = "
339  << fReactantsMV.size() << G4endl;
340 
341  G4int nbPrintable = fReactantsMV.size() * fReactantsMV.size();
342 
343  G4String* outputReaction = new G4String[nbPrintable];
344  G4String* outputReactionRate = new G4String[nbPrintable];
345  G4String* outputRange = new G4String[nbPrintable];
346  G4int n = 0;
347 
348  for (itReactives = fReactantsMV.begin(); itReactives != fReactantsMV.end();
349  itReactives++)
350  {
351  Reactant* moleculeA = (Reactant*)itReactives->first;
352  const vector<Reactant*>* reactivesVector = CanReactWith(moleculeA);
353 
354  if (pReactionModel) pReactionModel->InitialiseToPrint(moleculeA);
355 
356  G4int nbReactants = fReactantsMV[itReactives->first].size();
357 
358  for (G4int iReact = 0; iReact < nbReactants; iReact++)
359  {
360  auto moleculeB = (Reactant*)(*reactivesVector)[iReact];
361 
362  Data* reactionData = fReactionData[moleculeA][moleculeB];
363 
364  //-----------------------------------------------------------
365  // Name of the reaction
366  if (!alreadyPrint[moleculeA][moleculeB])
367  {
368  outputReaction[n] = moleculeA->GetName() + " + " + moleculeB->GetName();
369 
370  G4int nbProducts = reactionData->GetNbProducts();
371 
372  if (nbProducts)
373  {
374  outputReaction[n] += " -> " + reactionData->GetProduct(0)->GetName();
375 
376  for (G4int j = 1; j < nbProducts; j++)
377  {
378  outputReaction[n] += " + " + reactionData->GetProduct(j)->GetName();
379  }
380  }
381  else
382  {
383  outputReaction[n] += " -> No product";
384  }
385 
386  //-----------------------------------------------------------
387  // Interaction Rate
388  outputReactionRate[n] = G4UIcommand::ConvertToString(
389  reactionData->GetObservedReactionRateConstant() / (1e-3 * m3 / (mole * s)));
390 
391  //-----------------------------------------------------------
392  // Calculation of the Interaction Range
393  G4double interactionRange = -1;
394  if (pReactionModel) interactionRange =
395  pReactionModel->GetReactionRadius(iReact);
396 
397  if (interactionRange != -1)
398  {
399  outputRange[n] = G4UIcommand::ConvertToString(
400  interactionRange / nanometer);
401  }
402  else
403  {
404  outputRange[n] = "";
405  }
406 
407  alreadyPrint[moleculeB][moleculeA] = TRUE;
408  n++;
409  }
410  }
411  }
412  // G4cout<<"Number of possible reactions: "<< n << G4endl;
413 
415  // Tableau dynamique en fonction du nombre de caractere maximal dans
416  // chaque colonne
418 
419  G4int maxlengthOutputReaction = -1;
420  G4int maxlengthOutputReactionRate = -1;
421 
422  for (G4int i = 0; i < n; i++)
423  {
424  if (maxlengthOutputReaction < (G4int)outputReaction[i].length())
425  {
426  maxlengthOutputReaction = outputReaction[i].length();
427  }
428  if (maxlengthOutputReactionRate < (G4int)outputReactionRate[i].length())
429  {
430  maxlengthOutputReactionRate = outputReactionRate[i].length();
431  }
432  }
433 
434  maxlengthOutputReaction += 2;
435  maxlengthOutputReactionRate += 2;
436 
437  if (maxlengthOutputReaction < 10) maxlengthOutputReaction = 10;
438  if (maxlengthOutputReactionRate < 30) maxlengthOutputReactionRate = 30;
439 
440  G4String* title;
441 
442  if (pReactionModel) title = new G4String[3];
443  else title = new G4String[2];
444 
445  title[0] = "Reaction";
446  title[1] = "Reaction Rate [dm3/(mol*s)]";
447 
448  if (pReactionModel) title[2] =
449  "Interaction Range for chosen reaction model [nm]";
450 
451  G4cout << setfill(' ') << setw(maxlengthOutputReaction) << left << title[0]
452  << setw(maxlengthOutputReactionRate) << left << title[1];
453 
454  if (pReactionModel) G4cout << setw(2) << left << title[2];
455 
456  G4cout << G4endl;
457 
458  G4cout.fill('-');
459  if (pReactionModel) G4cout.width(
460  maxlengthOutputReaction + 2 + maxlengthOutputReactionRate + 2
461  + (G4int)title[2].length());
462  else G4cout.width(maxlengthOutputReaction + 2 + maxlengthOutputReactionRate);
463  G4cout << "-" << G4endl;
464  G4cout.fill(' ');
465 
466  for (G4int i = 0; i < n; i++)
467  {
468  G4cout << setw(maxlengthOutputReaction) << left << outputReaction[i]
469  << setw(maxlengthOutputReactionRate) << left
470  << outputReactionRate[i];
471 
472  if (pReactionModel) G4cout << setw(2) << left << outputRange[i];
473 
474  G4cout << G4endl;
475 
476  G4cout.fill('-');
477  if (pReactionModel) G4cout.width(
478  maxlengthOutputReaction + 2 + maxlengthOutputReactionRate + 2
479  + (G4int)title[2].length());
480  else G4cout.width(
481  maxlengthOutputReaction + 2 + maxlengthOutputReactionRate);
482  G4cout << "-" << G4endl;
483  G4cout.fill(' ');
484  }
485 
486  delete[] title;
487  delete[] outputReaction;
488  delete[] outputReactionRate;
489  delete[] outputRange;
490 }
491 
492 //______________________________________________________________________________
493 // Get/Set methods
494 
497  Reactant* pReactant2) const
498 {
499  if (fReactionData.empty())
500  {
501  G4String errMsg = "No reaction table was implemented";
502  G4Exception("G4MolecularInteractionTable::GetReactionData", "",
503  FatalErrorInArgument, errMsg);
504  }
505 
506  auto it1 = fReactionData.find(pReactant1);
507 
508  if (it1 == fReactionData.end())
509  {
510  G4String errMsg =
511  "No reaction table was implemented for this molecule Definition : " + pReactant1
512  ->GetName();
513  G4Exception("G4MolecularInteractionTable::GetReactionData", "",
514  FatalErrorInArgument, errMsg);
515  }
516 
517  ReactionDataMap::mapped_type::const_iterator it2 = it1->second.find(pReactant2);
518 
519  if (it2 == it1->second.end())
520  {
521  G4cout << "Name : " << pReactant2->GetName() << G4endl;
522  G4String errMsg = "No reaction table was implemented for this molecule : "
523  + pReactant2->GetName();
524  G4Exception("G4MolecularInteractionTable::GetReactionData", "", FatalErrorInArgument, errMsg);
525  }
526 
527  return (it2->second);
528 }
529 
531 {
532  return fReactionData;
533 }
534 
536 {
537  DataList dataList;
538 
539  for (const auto& pData : fVectorOfReactionData)
540  {
541  dataList.emplace_back(pData.get());
542  }
543 
544  return dataList;
545 }
546 
547 //______________________________________________________________________________
548 
551 {
552  if (fReactantsMV.empty())
553  {
554  G4String errMsg = "No reaction table was implemented";
555  G4Exception("G4MolecularInteractionTable::CanReactWith", "",
556  FatalErrorInArgument, errMsg);
557  return 0;
558  }
559 
560  auto itReactivesMap = fReactantsMV.find(pMolecule);
561 
562  if (itReactivesMap == fReactantsMV.end())
563  {
564 #ifdef G4VERBOSE
565  if (fVerbose)
566  {
567  G4String errMsg = "No reaction table was implemented for this molecule : "
568  + pMolecule->GetName();
569  // G4Exception("G4MolecularInteractionTable::CanReactWith","",FatalErrorInArgument, errMsg);
570  G4cout << "--- G4MolecularInteractionTable::GetReactionData ---" << G4endl;
571  G4cout << errMsg << G4endl;
572  }
573 #endif
574  return nullptr;
575  }
576 
577  if (fVerbose)
578  {
579  G4cout << " G4MolecularInteractionTable::CanReactWith :" << G4endl;
580  G4cout << "You are checking reactants for : " << pMolecule->GetName() << G4endl;
581  G4cout << " the number of reactants is : " << itReactivesMap->second.size() << G4endl;
582 
583  auto itProductsVector = itReactivesMap->second.cbegin();
584 
585  for (; itProductsVector != itReactivesMap->second.end(); itProductsVector++)
586  {
587  G4cout << (*itProductsVector)->GetName() << G4endl;
588  }
589  }
590  return &(itReactivesMap->second);
591 }
592 
593 //______________________________________________________________________________
594 
597 {
598  if (fReactionData.empty())
599  {
600  G4String errMsg = "No reaction table was implemented";
601  G4Exception("G4MolecularInteractionTable::CanInteractWith", "",
602  FatalErrorInArgument, errMsg);
603  }
604 
605  ReactionDataMap::const_iterator itReactivesMap = fReactionData.find(molecule);
606 
607  if (itReactivesMap == fReactionData.end())
608  {
609  return nullptr;
610  }
611 
612  if (fVerbose)
613  {
614  G4cout << " G4MolecularInteractionTable::CanReactWith :" << G4endl;
615  G4cout << "You are checking reactants for : " << molecule->GetName() << G4endl;
616  G4cout << " the number of reactants is : " << itReactivesMap->second.size() << G4endl;
617 
618  SpecificDataList::const_iterator itProductsVector = itReactivesMap->second.begin();
619 
620  for (; itProductsVector != itReactivesMap->second.end(); itProductsVector++)
621  {
622  G4cout << itProductsVector->first->GetName() << G4endl;
623  }
624  }
625  return &(itReactivesMap->second);
626 }
627 
628 //______________________________________________________________________________
629 
632 {
633  if (fReactionDataMV.empty())
634  {
635  G4String errMsg = "No reaction table was implemented";
636  G4Exception("G4MolecularInteractionTable::CanInteractWith", "",
637  FatalErrorInArgument, errMsg);
638  }
639  auto it = fReactionDataMV.find(molecule);
640 
641  if (it == fReactionDataMV.end())
642  {
643  G4String errMsg = "No reaction table was implemented for this molecule Definition : "
644  + molecule->GetName();
645  G4Exception("G4MolecularInteractionTable::GetReactionData", "", FatalErrorInArgument, errMsg);
646  }
647 
648  return &(it->second);
649 }
650 
651 //______________________________________________________________________________
652 
654  const G4String& mol2) const
655 {
656  const auto pConf1 = G4MoleculeTable::GetMoleculeTable()->GetConfiguration(mol1);
657  const auto pConf2 = G4MoleculeTable::GetMoleculeTable()->GetConfiguration(mol2);
658  return GetReactionData(pConf1, pConf2);
659 }
660 
661 //______________________________________________________________________________
662 
663 void
665 {
666  fRateParam = std::bind(PolynomialParam, std::placeholders::_1, P);
667 }
668 
669 //______________________________________________________________________________
670 
672  double E_R)
673 {
674  std::vector<double> P = { A0, E_R };
675  fRateParam = std::bind(ArrehniusParam, std::placeholders::_1, P);
676 }
677 
678 //______________________________________________________________________________
679 
681  double rateCste)
682 {
684  std::placeholders::_1,
685  temperature_K,
686  rateCste);
687 }
688 
689 //______________________________________________________________________________
690 
692 {
693  for (const auto& pData : fVectorOfReactionData)
694  {
695  const_cast<G4DNAMolecularReactionData*>(pData.get())->ScaleForNewTemperature(temp_K);
696  }
697 }
698 
699 //______________________________________________________________________________
700 
702 {
703  if (fRateParam)
704  {
706  }
707 }
708 
709 //______________________________________________________________________________
710 
713 {
714  for (auto& pData : fVectorOfReactionData)
715  {
716  if (pData->GetReactionID() == reactionID)
717  {
718  return pData.get();
719  }
720  }
721  return nullptr;
722 }
723 
725 {
726  return fVectorOfReactionData.size();
727 }