| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // Author: Mathieu Karamitros (kara (AT) cenbg
28 //
29 // WARNING : This class is released as a proto
30 // It might strongly evolve or even disapear i
31 //
32 // History:
33 // -----------
34 // 10 Oct 2011 M.Karamitros created
35 //
36 // -------------------------------------------
37
38 #include <iomanip>
39
40 #include "G4DNAMolecularReactionTable.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4UIcommand.hh"
44 #include "G4VDNAReactionModel.hh"
45 #include "G4MoleculeHandleManager.hh"
46 #include "G4MoleculeTable.hh"
47 #include "G4MolecularConfiguration.hh"
48 #include "G4ReactionTableMessenger.hh"
49 #include "G4IosFlagsSaver.hh"
50 #include "G4Exp.hh"
51
52 using namespace std;
53
54 G4DNAMolecularReactionTable* G4DNAMolecularRea
55
56 G4DNAMolecularReactionData::G4DNAMolecularReac
57 : fpReactant1(nullptr)
58 , fpReactant2(nullptr)
59 , fObservedReactionRate(0.)
60 , fActivationRate(0.)
61 , fDiffusionRate(0.)
62 , fOnsagerRadius(0.)
63 , fReactionRadius(0.)
64 , fEffectiveReactionRadius(0.)
65 , fProbability(0.)
66 , fType(0)
67 , fReactionID(0)
68 {
69 }
70
71 G4DNAMolecularReactionData::G4DNAMolecularReac
72
73
74 : fpReactant1(pReactant1)
75 , fpReactant2(pReactant2)
76 , fObservedReactionRate(reactionRate)
77 , fActivationRate(0.)
78 , fDiffusionRate(0.)
79 , fOnsagerRadius(0.)
80 , fReactionRadius(0.)
81 , fEffectiveReactionRadius(0.)
82 , fProbability(0.)
83 , fType(0)
84 , fReactionID(0)
85 {
86 ComputeEffectiveRadius();
87 }
88
89 G4DNAMolecularReactionData::G4DNAMolecularReac
90
91
92 : fpReactant1(nullptr)
93 , fpReactant2(nullptr)
94 , fObservedReactionRate(reactionRate)
95 , fActivationRate(0.)
96 , fDiffusionRate(0.)
97 , fOnsagerRadius(0.)
98 , fReactionRadius(0.)
99 , fEffectiveReactionRadius(0.)
100 , fProbability(0.)
101 , fType(0)
102 , fReactionID(0)
103 {
104 SetReactant1(reactant1);
105 SetReactant2(reactant2);
106 ComputeEffectiveRadius();
107 }
108
109 G4DNAMolecularReactionData::~G4DNAMolecularRea
110 {
111 fProducts.clear();
112 }
113
114 void G4DNAMolecularReactionData::ComputeEffect
115 {
116 G4double sumDiffCoeff = 0.;
117
118 if (fpReactant1 == fpReactant2)
119 {
120 sumDiffCoeff = fpReactant1->GetDiffusi
121 fEffectiveReactionRadius = fObservedRe
122 }
123 else
124 {
125 sumDiffCoeff = fpReactant1->GetDiffusi
126 + fpReactant2->GetDiffusi
127 fEffectiveReactionRadius = fObservedRe
128 }
129
130 fReactionID = 0;
131 fReactionRadius = fEffectiveReactionRadius
132 fOnsagerRadius = (fpReactant1->GetCharge()
133 fProbability = 1;
134
135 }
136
137 int G4DNAMolecularReactionData::GetReactionID(
138 {
139 return fReactionID;
140 }
141
142 void G4DNAMolecularReactionData::SetReactionID
143 {
144 fReactionID = ID;
145 }
146
147 void G4DNAMolecularReactionData::SetReactant1(
148 {
149 fpReactant1 = pReactive;
150 }
151
152 void G4DNAMolecularReactionData::SetReactant2(
153 {
154 fpReactant2 = pReactive;
155 }
156
157 void G4DNAMolecularReactionData::SetReactants(
158
159 {
160 fpReactant1 = pReactant1;
161 fpReactant2 = pReactant2;
162 }
163
164 void G4DNAMolecularReactionData::AddProduct(Re
165 {
166 fProducts.push_back(pMolecule);
167 }
168
169 G4int G4DNAMolecularReactionData::GetNbProduct
170 {
171 return (G4int)fProducts.size();
172 }
173
174 G4DNAMolecularReactionData::Reactant* G4DNAMol
175 {
176 return fProducts[i];
177 }
178
179 const G4DNAMolecularReactionData::ReactionProd
180 {
181 return &fProducts;
182 }
183
184 void G4DNAMolecularReactionData::RemoveProduct
185 {
186 fProducts.clear();
187 }
188
189 void G4DNAMolecularReactionData::SetReactant1(
190 {
191 fpReactant1 = G4MoleculeTable::Instance()-
192 }
193 void G4DNAMolecularReactionData::SetReactant2(
194 {
195 fpReactant2 = G4MoleculeTable::Instance()-
196 }
197 void G4DNAMolecularReactionData::SetReactants(
198
199 {
200 fpReactant1 = G4MoleculeTable::Instance()-
201 fpReactant2 = G4MoleculeTable::Instance()-
202 }
203
204 G4DNAMolecularReactionData::ReactantPair G4DNA
205 {
206 return std::make_pair(fpReactant1, fpReact
207 }
208
209 G4DNAMolecularReactionData::Reactant* G4DNAMol
210 {
211 return fpReactant1;
212 }
213
214 G4DNAMolecularReactionData::Reactant* G4DNAMol
215 {
216 return fpReactant2;
217 }
218
219 void G4DNAMolecularReactionData::SetObservedRe
220 {
221 fObservedReactionRate = rate;
222 }
223
224 G4double G4DNAMolecularReactionData::GetObserv
225 {
226 return fObservedReactionRate;
227 }
228
229 G4double G4DNAMolecularReactionData::GetActiva
230 {
231 return fActivationRate;
232 }
233
234 G4double G4DNAMolecularReactionData::GetDiffus
235 {
236 return fDiffusionRate;
237 }
238
239 void G4DNAMolecularReactionData::SetReactionRa
240 {
241 fReactionRadius = radius;
242 fEffectiveReactionRadius = -fOnsagerRadius
243 }
244
245 G4double G4DNAMolecularReactionData::GetReacti
246 {
247 return fReactionRadius;
248 }
249
250 void G4DNAMolecularReactionData::SetEffectiveR
251 {
252 fEffectiveReactionRadius = radius;
253 }
254
255 G4double G4DNAMolecularReactionData::GetEffect
256 {
257 return fEffectiveReactionRadius;
258 }
259
260 G4double G4DNAMolecularReactionData::GetOnsage
261 {
262 return fOnsagerRadius;
263 }
264
265 G4double G4DNAMolecularReactionData::GetProbab
266 {
267 return fProbability;
268 }
269
270 void G4DNAMolecularReactionData::SetProbabilit
271 {
272 fProbability = prob;
273 }
274
275 void G4DNAMolecularReactionData::SetReactionTy
276 {
277 G4double sumDiffCoeff = 0.;
278
279 if(type == 1)
280 {
281
282 sumDiffCoeff = fpReactant1->GetDiffusi
283 fpReactant2->GetDiffusi
284
285 fReactionRadius = fpReactant1->GetVanD
286 fpReactant2->GetVanD
287
288 G4double Rs = 0.29 * nm;
289 if(fOnsagerRadius == 0) // Type II
290 {
291 fEffectiveReactionRadius = fReacti
292 fDiffusionRate = 4 * pi * sumDiffC
293 if (fpReactant1 == fpReactant2) fD
294 fActivationRate = fDiffusionRate *
295 fProbability = Rs / (Rs + (fDiffu
296
297 }else{ // Type IV
298 fEffectiveReactionRadius = -fOnsag
299 fDiffusionRate = 4 * pi * sumDiffC
300 if (fpReactant1 == fpReactant2) fD
301
302 fActivationRate = fDiffusionRate *
303 fProbability = Rs / (Rs + (fDiffus
304 }
305 }
306
307 fType = type;
308 }
309
310 G4int G4DNAMolecularReactionData::GetReactionT
311 {
312 return fType;
313 }
314
315 void G4DNAMolecularReactionData::AddProduct(co
316 {
317 fProducts.push_back(G4MoleculeTable::Insta
318 }
319
320 double G4DNAMolecularReactionData::PolynomialP
321 {
322 double inv_temp = 1. / temp_K;
323
324 return pow(10,
325 P[0] + P[1] * inv_temp + P[2] *
326 + P[3] * pow(inv_temp, 3) + P[4
327 * (1e-3 * CLHEP::m3 / (CLHEP::mole * C
328 }
329
330 double G4DNAMolecularReactionData::ArrehniusPa
331 {
332 return P[0] * G4Exp(P[1] / temp_K)*
333 (1e-3 * CLHEP::m3 / (CLHEP::mole * CLH
334 }
335
336 double G4DNAMolecularReactionData::ScaledParam
337
338
339 {
340 double D0 = G4MolecularConfiguration::Diff
341 double Df = G4MolecularConfiguration::Diff
342 return Df * rateCste_init / D0;
343 }
344
345 //============================================
346 // REACTION TABLE
347 //============================================
348
349 G4DNAMolecularReactionTable* G4DNAMolecularRea
350 {
351 if (fpInstance == nullptr)
352 {
353 fpInstance = new G4DNAMolecularReactio
354 }
355 return fpInstance;
356 }
357
358 //____________________________________________
359
360 G4DNAMolecularReactionTable* G4DNAMolecularRea
361 {
362 if (fpInstance == nullptr)
363 {
364 fpInstance = new G4DNAMolecularReactio
365 }
366 return fpInstance;
367 }
368
369 //____________________________________________
370
371 void G4DNAMolecularReactionTable::DeleteInstan
372 {
373
374
375 delete fpInstance;
376
377 fpInstance = nullptr;
378 }
379
380 //____________________________________________
381
382 G4DNAMolecularReactionTable::G4DNAMolecularRea
383 :
384 fpMessenger(new G4ReactionTableMessenger(
385 {
386 }
387
388 G4DNAMolecularReactionTable::~G4DNAMolecularRe
389
390 void G4DNAMolecularReactionTable::SetReaction(
391 {
392 const auto pReactant1 = pReactionData->Get
393 const auto pReactant2 = pReactionData->Get
394
395 fReactionData[pReactant1][pReactant2] = pR
396 fReactantsMV[pReactant1].push_back(pReacta
397 fReactionDataMV[pReactant1].push_back(pRea
398
399 if (pReactant1 != pReactant2)
400 {
401 fReactionData[pReactant2][pReactant1]
402 fReactantsMV[pReactant2].push_back(pRe
403 fReactionDataMV[pReactant2].push_back(
404 }
405
406 fVectorOfReactionData.emplace_back(pReacti
407 pReactionData->SetReactionID((G4int)fVecto
408 }
409
410 //____________________________________________
411
412 void G4DNAMolecularReactionTable::SetReaction(
413
414
415 {
416 auto reactionData = new G4DNAMolecularReac
417 SetReaction(reactionData);
418 }
419
420 //____________________________________________
421
422 void G4DNAMolecularReactionTable::PrintTable(G
423 {
424 // Print Reactions and Interaction radius
425
426 G4IosFlagsSaver iosfs(G4cout);
427
428 if ((pReactionModel != nullptr) && ((pReact
429 {
430 pReactionModel->SetReactionTable(this);
431 }
432
433 ReactivesMV::iterator itReactives;
434
435 std::map<Reactant*, std::map<Reactant*, G4
436
437 G4cout << "Number of chemical species invo
438 << fReactantsMV.size() << G4endl;
439
440 std::size_t nbPrintable = fReactantsMV.siz
441
442 auto outputReaction = new G4String[nbPrin
443 auto outputReactionRate = new G4String[nb
444 auto outputRange = new G4String[nbPrintab
445 G4int n = 0;
446
447 for (itReactives = fReactantsMV.begin(); i
448 ++itReactives)
449 {
450 auto moleculeA = (Reactant*)itReactiv
451 const vector<Reactant*>* reactivesVect
452
453 if (pReactionModel != nullptr) pReacti
454
455 auto nbReactants = (G4int)fReactantsM
456
457 for (G4int iReact = 0; iReact < nbReac
458 {
459 auto moleculeB = (Reactant*)(*reac
460
461 Data* reactionData = fReactionData
462
463 //--------------------------------
464 // Name of the reaction
465 if (!alreadyPrint[moleculeA][molec
466 {
467 outputReaction[n] = moleculeA-
468
469 G4int nbProducts = reactionDat
470
471 if (nbProducts != 0)
472 {
473 outputReaction[n] += " ->
474
475 for (G4int j = 1; j < nbPr
476 {
477 outputReaction[n] += "
478 }
479 }
480 else
481 {
482 outputReaction[n] += " ->
483 }
484
485 //----------------------------
486 // Interaction Rate
487 outputReactionRate[n] = G4UIco
488 reactionData->GetObservedR
489
490 //----------------------------
491 // Calculation of the Interact
492 G4double interactionRange = -1
493 if (pReactionModel != nullptr)
494 pReactionModel->GetReactio
495
496 if (interactionRange != -1)
497 {
498 outputRange[n] = G4UIcomma
499 interactionRange / nan
500 }
501 else
502 {
503 outputRange[n] = "";
504 }
505
506 alreadyPrint[moleculeB][molecu
507 n++;
508 }
509 }
510 }
511 // G4cout<<"Number of possible reactions:
512
513 //////////////////////////////////////////
514 // Tableau dynamique en fonction du nombre
515 // chaque colonne
516 //////////////////////////////////////////
517
518 G4int maxlengthOutputReaction = -1;
519 G4int maxlengthOutputReactionRate = -1;
520
521 for (G4int i = 0; i < n; ++i)
522 {
523 if (maxlengthOutputReaction < (G4int)o
524 {
525 maxlengthOutputReaction = (G4int)o
526 }
527 if (maxlengthOutputReactionRate < (G4i
528 {
529 maxlengthOutputReactionRate = (G4i
530 }
531 }
532
533 maxlengthOutputReaction += 2;
534 maxlengthOutputReactionRate += 2;
535
536 if (maxlengthOutputReaction < 10) maxlengt
537 if (maxlengthOutputReactionRate < 30) maxl
538
539 G4String* title;
540
541 if (pReactionModel != nullptr) title = new
542 else title = new G4String[2];
543
544 title[0] = "Reaction";
545 title[1] = "Reaction Rate [dm3/(mol*s)]";
546
547 if (pReactionModel != nullptr) title[2] =
548 "Interaction Range for chosen reaction
549
550 G4cout << setfill(' ') << setw(maxlengthOu
551 << setw(maxlengthOutputReactionRate) <
552
553 if (pReactionModel != nullptr) G4cout << s
554
555 G4cout << G4endl;
556
557 G4cout.fill('-');
558 if (pReactionModel != nullptr) G4cout.widt
559 maxlengthOutputReaction + 2 + maxlengt
560 + (G4int)title[2].length());
561 else G4cout.width(maxlengthOutputReaction
562 G4cout << "-" << G4endl;
563 G4cout.fill(' ');
564
565 for (G4int i = 0; i < n; i++)
566 {
567 G4cout << setw(maxlengthOutputReaction
568 << setw(maxlengthOutputReaction
569 << outputReactionRate[i];
570
571 if (pReactionModel != nullptr) G4cout
572
573 G4cout << G4endl;
574
575 G4cout.fill('-');
576 if (pReactionModel != nullptr) G4cout.
577 maxlengthOutputReaction + 2 + maxl
578 + (G4int)title[2].length());
579 else G4cout.width(
580 maxlengthOutputReaction + 2 + maxl
581 G4cout << "-" << G4endl;
582 G4cout.fill(' ');
583 }
584
585 delete[] title;
586 delete[] outputReaction;
587 delete[] outputReactionRate;
588 delete[] outputRange;
589 }
590
591 //____________________________________________
592 // Get/Set methods
593
594 G4VDNAMolecularGeometry* G4DNAMolecularReactio
595 {
596 return fGeometry;
597 }
598
599 G4DNAMolecularReactionTable::Data*
600 G4DNAMolecularReactionTable::GetReactionData(R
601 R
602 {
603 if (fReactionData.empty())
604 {
605 G4String errMsg = "No reaction table w
606 G4Exception("G4MolecularInteractionTab
607 FatalErrorInArgument, errM
608 }
609
610 auto it1 = fReactionData.find(pReactant1);
611
612 if (it1 == fReactionData.end())
613 {
614 G4String errMsg =
615 "No reaction table was implemented
616 ->GetName();
617 G4Exception("G4MolecularInteractionTab
618 FatalErrorInArgument, errM
619 // Though the above is Fatal and will
620 return nullptr;
621 }
622
623 auto it2 = it1->second.find(pReactant2);
624
625 if (it2 == it1->second.end())
626 {
627 G4cout << "Name : " << pReactant2->Get
628 G4String errMsg = "No reaction table w
629 + pReactant2->GetName();
630 G4Exception("G4MolecularInteractionTab
631 }
632
633 return (it2->second);
634 }
635
636 const G4DNAMolecularReactionTable::ReactionDat
637 {
638 return fReactionData;
639 }
640
641 G4DNAMolecularReactionTable::DataList G4DNAMol
642 {
643 DataList dataList;
644
645 for (const auto& pData : fVectorOfReaction
646 {
647 dataList.emplace_back(pData.get());
648 }
649
650 return dataList;
651 }
652
653 //____________________________________________
654
655 const G4DNAMolecularReactionTable::ReactantLis
656 G4DNAMolecularReactionTable::CanReactWith(Reac
657 {
658 if (fReactantsMV.empty())
659 {
660 G4String errMsg = "No reaction table w
661 G4Exception("G4MolecularInteractionTab
662 FatalErrorInArgument, errM
663 return nullptr;
664 }
665
666 auto itReactivesMap = fReactantsMV.find(pM
667
668 if (itReactivesMap == fReactantsMV.end())
669 {
670 #ifdef G4VERBOSE
671 if (fVerbose)
672 {
673 G4String errMsg = "No reaction tab
674 + pMolecule->GetName();
675 // G4Exception("G4Molecul
676 G4cout << "--- G4MolecularInteract
677 G4cout << errMsg << G4endl;
678 }
679 #endif
680 return nullptr;
681 }
682
683 if (fVerbose)
684 {
685 G4cout << " G4MolecularInteractionTabl
686 G4cout << "You are checking reactants
687 G4cout << " the number of reactants is
688
689 auto itProductsVector = itReactivesMap
690
691 for (; itProductsVector != itReactives
692 {
693 G4cout << (*itProductsVector)->Get
694 }
695 }
696 return &(itReactivesMap->second);
697 }
698
699 //____________________________________________
700
701 const G4DNAMolecularReactionTable::SpecificDat
702 G4DNAMolecularReactionTable::GetReativesNData(
703 {
704 if (fReactionData.empty())
705 {
706 G4String errMsg = "No reaction table w
707 G4Exception("G4MolecularInteractionTab
708 FatalErrorInArgument, errM
709 }
710
711 auto itReactivesMap = fReactionData.find(m
712
713 if (itReactivesMap == fReactionData.end())
714 {
715 return nullptr;
716 }
717
718 if (fVerbose)
719 {
720 G4cout << " G4MolecularInteractionTabl
721 G4cout << "You are checking reactants
722 G4cout << " the number of reactants is
723
724 auto itProductsVector = itReactivesMap
725
726 for (; itProductsVector != itReactives
727 {
728 G4cout << itProductsVector->first-
729 }
730 }
731 return &(itReactivesMap->second);
732 }
733
734 //____________________________________________
735
736 const G4DNAMolecularReactionTable::DataList*
737 G4DNAMolecularReactionTable::GetReactionData(c
738 {
739 if (fReactionDataMV.empty())
740 {
741 G4String errMsg = "No reaction table w
742 G4Exception("G4MolecularInteractionTab
743 FatalErrorInArgument, errM
744 }
745 auto it = fReactionDataMV.find(molecule);
746
747 if (it == fReactionDataMV.end())
748 {
749 G4String errMsg = "No reaction table w
750 + molecule->GetName();
751 G4Exception("G4MolecularInteractionTab
752 // Though the above is Fatal and will
753 return nullptr;
754 }
755
756 return &(it->second);
757 }
758
759 //____________________________________________
760
761 G4DNAMolecularReactionTable::Data* G4DNAMolecu
762
763 {
764 const auto pConf1 = G4MoleculeTable::GetMo
765 const auto pConf2 = G4MoleculeTable::GetMo
766 return GetReactionData(pConf1, pConf2);
767 }
768
769 //____________________________________________
770
771 void
772 G4DNAMolecularReactionData::SetPolynomialParam
773 {
774 fRateParam = std::bind(PolynomialParam, st
775 }
776
777 //____________________________________________
778
779 void G4DNAMolecularReactionData::SetArrehniusP
780
781 {
782 std::vector<double> P = { A0, E_R };
783 fRateParam = std::bind(ArrehniusParam, std
784 }
785
786 //____________________________________________
787
788 void G4DNAMolecularReactionData::SetScaledPara
789
790 {
791 fRateParam = std::bind(ScaledParameterizat
792 std::placeholders::
793 temperature_K,
794 rateCste);
795 }
796
797 //____________________________________________
798
799 void G4DNAMolecularReactionTable::ScaleReactio
800 {
801 for (const auto& pData : fVectorOfReaction
802 {
803 const_cast<G4DNAMolecularReactionData*
804 }
805 }
806
807 //____________________________________________
808
809 void G4DNAMolecularReactionData::ScaleForNewTe
810 {
811 if (fRateParam)
812 {
813 SetObservedReactionRateConstant(fRateP
814 }
815 }
816
817 //____________________________________________
818
819 G4DNAMolecularReactionTable::Data*
820 G4DNAMolecularReactionTable::GetReaction(int r
821 {
822 for (auto& pData : fVectorOfReactionData)
823 {
824 if (pData->GetReactionID() == reaction
825 {
826 return pData.get();
827 }
828 }
829 return nullptr;
830 }
831
832 size_t G4DNAMolecularReactionTable::GetNReacti
833 {
834 return fVectorOfReactionData.size();
835 }
836
837 void G4DNAMolecularReactionTable::Reset()
838 {
839 fReactionData.clear();
840 fReactantsMV.clear();
841 fReactionDataMV.clear();
842 fVectorOfReactionData.clear();
843 }
844