| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4DNAMolecularReactionTable.cc 65022 2012-11-12 16:43:12Z gcosmo $
26 // 27 //
27 // Author: Mathieu Karamitros (kara (AT) cenbg 28 // Author: Mathieu Karamitros (kara (AT) cenbg . in2p3 . fr)
28 // 29 //
29 // WARNING : This class is released as a proto 30 // WARNING : This class is released as a prototype.
30 // It might strongly evolve or even disapear i 31 // It might strongly evolve or even disapear in the next releases.
31 // 32 //
32 // History: 33 // History:
33 // ----------- 34 // -----------
34 // 10 Oct 2011 M.Karamitros created 35 // 10 Oct 2011 M.Karamitros created
35 // 36 //
36 // ------------------------------------------- 37 // -------------------------------------------------------------------
37 38
38 #include <iomanip> 39 #include <iomanip>
39 40
40 #include "G4DNAMolecularReactionTable.hh" 41 #include "G4DNAMolecularReactionTable.hh"
41 #include "G4PhysicalConstants.hh" 42 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh" 43 #include "G4SystemOfUnits.hh"
43 #include "G4UIcommand.hh" 44 #include "G4UIcommand.hh"
44 #include "G4VDNAReactionModel.hh" 45 #include "G4VDNAReactionModel.hh"
45 #include "G4MoleculeHandleManager.hh" 46 #include "G4MoleculeHandleManager.hh"
46 #include "G4MoleculeTable.hh" <<
47 #include "G4MolecularConfiguration.hh" <<
48 #include "G4ReactionTableMessenger.hh" <<
49 #include "G4IosFlagsSaver.hh" <<
50 #include "G4Exp.hh" <<
51 47
52 using namespace std; 48 using namespace std;
53 49
54 G4DNAMolecularReactionTable* G4DNAMolecularRea << 50 G4DNAMolecularReactionTable* G4DNAMolecularReactionTable::fInstance(0);
55 51
56 G4DNAMolecularReactionData::G4DNAMolecularReac << 52 G4DNAMolecularReactionData::G4DNAMolecularReactionData():
57 : fpReactant1(nullptr) << 53 fReactive1(),fReactive2(),
58 , fpReactant2(nullptr) << 54 fReactionRate(0.),fReducedReactionRadius(0.),
59 , fObservedReactionRate(0.) << 55 fProducts(0)
60 , fActivationRate(0.) << 56 {;}
61 , fDiffusionRate(0.) <<
62 , fOnsagerRadius(0.) <<
63 , fReactionRadius(0.) <<
64 , fEffectiveReactionRadius(0.) <<
65 , fProbability(0.) <<
66 , fType(0) <<
67 , fReactionID(0) <<
68 { <<
69 } <<
70 57
71 G4DNAMolecularReactionData::G4DNAMolecularReac 58 G4DNAMolecularReactionData::G4DNAMolecularReactionData(G4double reactionRate,
72 << 59 const G4Molecule* reactive1,
73 << 60 const G4Molecule* reactive2):fProducts(0)
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 { 61 {
86 ComputeEffectiveRadius(); << 62 fReactionRate = reactionRate;
87 } << 63 SetReactive1(reactive1);
>> 64 SetReactive2(reactive2);
88 65
89 G4DNAMolecularReactionData::G4DNAMolecularReac << 66 G4double sumDiffCoeff(0.);
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 67
109 G4DNAMolecularReactionData::~G4DNAMolecularRea << 68 if(*reactive1 == *reactive2)
110 { <<
111 fProducts.clear(); <<
112 } <<
113 <<
114 void G4DNAMolecularReactionData::ComputeEffect <<
115 { <<
116 G4double sumDiffCoeff = 0.; <<
117 <<
118 if (fpReactant1 == fpReactant2) <<
119 { 69 {
120 sumDiffCoeff = fpReactant1->GetDiffusi << 70 sumDiffCoeff = reactive1->GetDiffusionCoefficient();
121 fEffectiveReactionRadius = fObservedRe << 71 fReducedReactionRadius = fReactionRate/(4*pi* reactive1->GetDiffusionCoefficient() * Avogadro);
122 } 72 }
123 else 73 else
124 { 74 {
125 sumDiffCoeff = fpReactant1->GetDiffusi << 75 sumDiffCoeff = reactive1->GetDiffusionCoefficient() + reactive2->GetDiffusionCoefficient();
126 + fpReactant2->GetDiffusi << 76 fReducedReactionRadius = fReactionRate/(4*pi* sumDiffCoeff * Avogadro);
127 fEffectiveReactionRadius = fObservedRe <<
128 } 77 }
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 } 78 }
178 79
179 const G4DNAMolecularReactionData::ReactionProd << 80 G4DNAMolecularReactionData::~G4DNAMolecularReactionData()
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 { 81 {
252 fEffectiveReactionRadius = radius; << 82 if(fProducts)
>> 83 {
>> 84 fProducts->clear() ;
>> 85 delete fProducts;
>> 86 fProducts = 0;
>> 87 }
253 } 88 }
254 89
255 G4double G4DNAMolecularReactionData::GetEffect << 90 void G4DNAMolecularReactionData::SetReactive1(const G4Molecule* reactive)
256 { 91 {
257 return fEffectiveReactionRadius; << 92 fReactive1 = G4MoleculeHandleManager::Instance()->GetMoleculeHandle(reactive);
258 } 93 }
259 << 94 void G4DNAMolecularReactionData::SetReactive2(const G4Molecule* reactive)
260 G4double G4DNAMolecularReactionData::GetOnsage <<
261 { 95 {
262 return fOnsagerRadius; << 96 fReactive2 = G4MoleculeHandleManager::Instance()->GetMoleculeHandle(reactive);
263 } 97 }
264 << 98 void G4DNAMolecularReactionData::SetReactive(const G4Molecule* reactive1,
265 G4double G4DNAMolecularReactionData::GetProbab << 99 const G4Molecule* reactive2)
266 { 100 {
267 return fProbability; << 101 fReactive1 = G4MoleculeHandleManager::Instance()->GetMoleculeHandle(reactive1);
>> 102 fReactive2 = G4MoleculeHandleManager::Instance()->GetMoleculeHandle(reactive2);
268 } 103 }
269 104
270 void G4DNAMolecularReactionData::SetProbabilit << 105 void G4DNAMolecularReactionData::AddProduct(const G4Molecule* molecule)
271 { 106 {
272 fProbability = prob; << 107 if(!fProducts) fProducts = new std::vector<G4MoleculeHandle>();
>> 108 fProducts->push_back(G4MoleculeHandleManager::Instance()->GetMoleculeHandle(molecule));
273 } 109 }
274 << 110 //_____________________________________________________________________________________
275 void G4DNAMolecularReactionData::SetReactionTy << 111 G4DNAMolecularReactionTable* G4DNAMolecularReactionTable::GetReactionTable()
276 { 112 {
277 G4double sumDiffCoeff = 0.; << 113 if(!fInstance)
278 <<
279 if(type == 1) <<
280 { 114 {
281 << 115 fInstance = new G4DNAMolecularReactionTable();
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 } 116 }
306 << 117 return fInstance;
307 fType = type; <<
308 } <<
309 <<
310 G4int G4DNAMolecularReactionData::GetReactionT <<
311 { <<
312 return fType; <<
313 } 118 }
314 119
315 void G4DNAMolecularReactionData::AddProduct(co << 120 void G4DNAMolecularReactionTable::DeleteInstance()
316 { <<
317 fProducts.push_back(G4MoleculeTable::Insta <<
318 } <<
319 <<
320 double G4DNAMolecularReactionData::PolynomialP <<
321 { 121 {
322 double inv_temp = 1. / temp_K; << 122 // DEBUG
323 << 123 // G4cout << "G4MolecularReactionTable::DeleteInstance" << G4endl;
324 return pow(10, << 124 if(fInstance)
325 P[0] + P[1] * inv_temp + P[2] * << 125 delete fInstance;
326 + P[3] * pow(inv_temp, 3) + P[4 << 126 fInstance = 0;
327 * (1e-3 * CLHEP::m3 / (CLHEP::mole * C <<
328 } 127 }
329 << 128 //_____________________________________________________________________________________
330 double G4DNAMolecularReactionData::ArrehniusPa << 129 G4DNAMolecularReactionTable::G4DNAMolecularReactionTable() : G4ITReactionTable(),
>> 130 fMoleculeHandleManager(G4MoleculeHandleManager::Instance())
331 { 131 {
332 return P[0] * G4Exp(P[1] / temp_K)* << 132 // G4cout << "G4DNAMolecularReactionTable::G4DNAMolecularReactionTable()" << G4endl;
333 (1e-3 * CLHEP::m3 / (CLHEP::mole * CLH << 133 fVerbose = false;
>> 134 return;
334 } 135 }
335 << 136 //_____________________________________________________________________________________
336 double G4DNAMolecularReactionData::ScaledParam << 137 G4DNAMolecularReactionTable::~G4DNAMolecularReactionTable()
337 <<
338 <<
339 { 138 {
340 double D0 = G4MolecularConfiguration::Diff << 139 // DEBUG
341 double Df = G4MolecularConfiguration::Diff << 140 // G4cout << "G4MolecularReactionTable::~G4MolecularReactionTable" << G4endl;
342 return Df * rateCste_init / D0; << 141 ReactionDataMap::iterator it1 = fReactionData.begin();
343 } <<
344 142
345 //============================================ << 143 std::map<const G4Molecule*,
346 // REACTION TABLE << 144 const G4DNAMolecularReactionData*,
347 //============================================ << 145 compMoleculeP>::iterator it2;
348 146
349 G4DNAMolecularReactionTable* G4DNAMolecularRea << 147 for(;it1!=fReactionData.end();it1++)
350 { <<
351 if (fpInstance == nullptr) <<
352 { 148 {
353 fpInstance = new G4DNAMolecularReactio << 149 for(it2 = it1->second.begin();it2 != it1->second.end();it2++)
354 } << 150 {
355 return fpInstance; << 151 const G4DNAMolecularReactionData* reactionData = it2->second;
356 } << 152 if(reactionData)
>> 153 {
>> 154 const G4Molecule* reactive1 = reactionData->GetReactive1();
>> 155 const G4Molecule* reactive2 = reactionData->GetReactive2();
357 156
358 //____________________________________________ << 157 fReactionData[reactive1][reactive2] = 0;
>> 158 fReactionData[reactive2][reactive1] = 0;
359 159
360 G4DNAMolecularReactionTable* G4DNAMolecularRea << 160 delete reactionData;
361 { << 161 }
362 if (fpInstance == nullptr) << 162 }
363 { <<
364 fpInstance = new G4DNAMolecularReactio <<
365 } 163 }
366 return fpInstance; <<
367 } <<
368 <<
369 //____________________________________________ <<
370 164
371 void G4DNAMolecularReactionTable::DeleteInstan << 165 fReactionDataMV.clear();
372 { << 166 fReactionData.clear();
373 << 167 fReactivesMV.clear();
374 <<
375 delete fpInstance; <<
376 <<
377 fpInstance = nullptr; <<
378 } 168 }
379 <<
380 //____________________________________________ 169 //_____________________________________________________________________________________
381 << 170 void G4DNAMolecularReactionTable::SetReaction(G4DNAMolecularReactionData* reactionData)
382 G4DNAMolecularReactionTable::G4DNAMolecularRea <<
383 : <<
384 fpMessenger(new G4ReactionTableMessenger( <<
385 { <<
386 } <<
387 <<
388 G4DNAMolecularReactionTable::~G4DNAMolecularRe <<
389 <<
390 void G4DNAMolecularReactionTable::SetReaction( <<
391 { 171 {
392 const auto pReactant1 = pReactionData->Get << 172 const G4Molecule* reactive1 = reactionData->GetReactive1() ;
393 const auto pReactant2 = pReactionData->Get << 173 const G4Molecule* reactive2 = reactionData->GetReactive2() ;
394 174
395 fReactionData[pReactant1][pReactant2] = pR << 175 fReactionData[reactive1][reactive2] = reactionData;
396 fReactantsMV[pReactant1].push_back(pReacta << 176 fReactivesMV[reactive1].push_back(reactive2);
397 fReactionDataMV[pReactant1].push_back(pRea << 177 fReactionDataMV[reactive1].push_back(reactionData);
398 178
399 if (pReactant1 != pReactant2) << 179 if(reactive1 != reactive2)
400 { 180 {
401 fReactionData[pReactant2][pReactant1] << 181 fReactionData[reactive2][reactive1] = reactionData;
402 fReactantsMV[pReactant2].push_back(pRe << 182 fReactivesMV[reactive2].push_back(reactive1);
403 fReactionDataMV[pReactant2].push_back( << 183 fReactionDataMV[reactive2].push_back(reactionData);
404 } 184 }
405 <<
406 fVectorOfReactionData.emplace_back(pReacti <<
407 pReactionData->SetReactionID((G4int)fVecto <<
408 } 185 }
409 <<
410 //____________________________________________ 186 //_____________________________________________________________________________________
411 <<
412 void G4DNAMolecularReactionTable::SetReaction( 187 void G4DNAMolecularReactionTable::SetReaction(G4double reactionRate,
413 << 188 const G4Molecule* reactive1,
414 << 189 const G4Molecule* reactive2)
415 { 190 {
416 auto reactionData = new G4DNAMolecularReac << 191 G4DNAMolecularReactionData* reactionData = new G4DNAMolecularReactionData(reactionRate, reactive1, reactive2);
417 SetReaction(reactionData); 192 SetReaction(reactionData);
418 } 193 }
419 <<
420 //____________________________________________ 194 //_____________________________________________________________________________________
421 <<
422 void G4DNAMolecularReactionTable::PrintTable(G 195 void G4DNAMolecularReactionTable::PrintTable(G4VDNAReactionModel* pReactionModel)
423 { 196 {
424 // Print Reactions and Interaction radius 197 // Print Reactions and Interaction radius for jump step = 3ps
425 198
426 G4IosFlagsSaver iosfs(G4cout); << 199 if(pReactionModel)
427 << 200 {
428 if ((pReactionModel != nullptr) && ((pReact << 201 if(!(pReactionModel->GetReactionTable()))
429 { << 202 pReactionModel -> SetReactionTable(this);
430 pReactionModel->SetReactionTable(this); << 203 }
431 } <<
432 204
433 ReactivesMV::iterator itReactives; 205 ReactivesMV::iterator itReactives;
434 206
435 std::map<Reactant*, std::map<Reactant*, G4 << 207 map<G4Molecule*,map<G4Molecule*, G4bool> > alreadyPrint;
436 208
437 G4cout << "Number of chemical species invo << 209 G4cout<<"Nombre particules intervenants dans les reactions = "<< fReactivesMV.size() <<G4endl;
438 << fReactantsMV.size() << G4endl; <<
439 210
440 std::size_t nbPrintable = fReactantsMV.siz << 211 G4int nbPrintable = fReactivesMV.size()*fReactivesMV.size();
441 212
442 auto outputReaction = new G4String[nbPrin << 213 G4String *outputReaction = new G4String[nbPrintable];
443 auto outputReactionRate = new G4String[nb << 214 G4String *outputReactionRate = new G4String[nbPrintable];
444 auto outputRange = new G4String[nbPrintab << 215 G4String *outputRange = new G4String[nbPrintable];
445 G4int n = 0; 216 G4int n = 0;
446 217
447 for (itReactives = fReactantsMV.begin(); i << 218 for(itReactives = fReactivesMV.begin() ; itReactives != fReactivesMV.end() ; itReactives++)
448 ++itReactives) <<
449 { 219 {
450 auto moleculeA = (Reactant*)itReactiv << 220 G4Molecule* moleculeA = (G4Molecule*) itReactives->first;
451 const vector<Reactant*>* reactivesVect << 221 const vector<const G4Molecule*>* reactivesVector = CanReactWith(moleculeA);
452 222
453 if (pReactionModel != nullptr) pReacti << 223 if(pReactionModel)
>> 224 pReactionModel -> InitialiseToPrint(moleculeA);
454 225
455 auto nbReactants = (G4int)fReactantsM << 226 G4int nbReactants = fReactivesMV[itReactives->first].size();
456 227
457 for (G4int iReact = 0; iReact < nbReac << 228 for(G4int iReact = 0 ; iReact < nbReactants ; iReact++)
458 { 229 {
459 auto moleculeB = (Reactant*)(*reac <<
460 230
461 Data* reactionData = fReactionData << 231 G4Molecule* moleculeB = (G4Molecule*) (*reactivesVector)[iReact];
>> 232
>> 233 const G4DNAMolecularReactionData* reactionData = fReactionData[moleculeA][moleculeB];
462 234
463 //-------------------------------- 235 //-----------------------------------------------------------
464 // Name of the reaction 236 // Name of the reaction
465 if (!alreadyPrint[moleculeA][molec << 237 if(!alreadyPrint[moleculeA][moleculeB])
466 { 238 {
467 outputReaction[n] = moleculeA- << 239 outputReaction[n]=
>> 240 moleculeA->GetName()
>> 241 +" + " +
>> 242 moleculeB->GetName();
468 243
469 G4int nbProducts = reactionDat 244 G4int nbProducts = reactionData->GetNbProducts();
470 245
471 if (nbProducts != 0) << 246 if(nbProducts)
472 { 247 {
473 outputReaction[n] += " -> << 248 outputReaction[n] += " -> "+ reactionData->GetProduct(0)->GetName();
474 249
475 for (G4int j = 1; j < nbPr << 250 for(G4int j = 1 ; j < nbProducts ; j++)
476 { 251 {
477 outputReaction[n] += " << 252 outputReaction[n]+=" + "+reactionData->GetProduct(j)->GetName();
478 } 253 }
479 } 254 }
480 else 255 else
481 { 256 {
482 outputReaction[n] += " -> << 257 outputReaction[n]+=" -> No product";
483 } 258 }
484 259
485 //---------------------------- 260 //-----------------------------------------------------------
486 // Interaction Rate 261 // Interaction Rate
487 outputReactionRate[n] = G4UIco << 262 outputReactionRate[n] = G4UIcommand::ConvertToString(reactionData->GetReactionRate()/(1e-3*m3/(mole*s)));
488 reactionData->GetObservedR <<
489 263
490 //---------------------------- 264 //-----------------------------------------------------------
491 // Calculation of the Interact 265 // Calculation of the Interaction Range
492 G4double interactionRange = -1 266 G4double interactionRange = -1;
493 if (pReactionModel != nullptr) << 267 if(pReactionModel)
494 pReactionModel->GetReactio << 268 interactionRange = pReactionModel->GetReactionRadius(iReact);
495 269
496 if (interactionRange != -1) << 270 if(interactionRange!=-1)
497 { 271 {
498 outputRange[n] = G4UIcomma << 272 outputRange[n] = G4UIcommand::ConvertToString(interactionRange/nanometer);
499 interactionRange / nan <<
500 } 273 }
501 else 274 else
502 { 275 {
503 outputRange[n] = ""; 276 outputRange[n] = "";
504 } 277 }
505 278
506 alreadyPrint[moleculeB][molecu 279 alreadyPrint[moleculeB][moleculeA] = TRUE;
507 n++; 280 n++;
508 } 281 }
509 } 282 }
510 } 283 }
511 // G4cout<<"Number of possible reactions: << 284 G4cout<<"Number of possible reactions: "<< n << G4endl;
512 285
513 ////////////////////////////////////////// 286 ////////////////////////////////////////////////////////////////////
514 // Tableau dynamique en fonction du nombre << 287 // Tableau dynamique en fonction du nombre de caractère maximal dans
515 // chaque colonne 288 // chaque colonne
516 ////////////////////////////////////////// 289 ////////////////////////////////////////////////////////////////////
517 290
518 G4int maxlengthOutputReaction = -1; 291 G4int maxlengthOutputReaction = -1;
519 G4int maxlengthOutputReactionRate = -1; 292 G4int maxlengthOutputReactionRate = -1;
520 293
521 for (G4int i = 0; i < n; ++i) << 294 for(G4int i = 0 ; i < n ; i++)
522 { 295 {
523 if (maxlengthOutputReaction < (G4int)o << 296 if(maxlengthOutputReaction < (G4int) outputReaction[i].length())
524 { 297 {
525 maxlengthOutputReaction = (G4int)o << 298 maxlengthOutputReaction = outputReaction[i].length();
526 } 299 }
527 if (maxlengthOutputReactionRate < (G4i << 300 if(maxlengthOutputReactionRate < (G4int)outputReactionRate[i].length())
528 { 301 {
529 maxlengthOutputReactionRate = (G4i << 302 maxlengthOutputReactionRate = outputReactionRate[i].length();
530 } 303 }
531 } 304 }
532 305
533 maxlengthOutputReaction += 2; << 306 maxlengthOutputReaction+=2;
534 maxlengthOutputReactionRate += 2; << 307 maxlengthOutputReactionRate+=2;
535 308
536 if (maxlengthOutputReaction < 10) maxlengt << 309 if(maxlengthOutputReaction<10) maxlengthOutputReaction = 10;
537 if (maxlengthOutputReactionRate < 30) maxl << 310 if(maxlengthOutputReactionRate<30) maxlengthOutputReactionRate = 30;
538 311
539 G4String* title; << 312 G4String title[3];
540 <<
541 if (pReactionModel != nullptr) title = new <<
542 else title = new G4String[2]; <<
543 313
544 title[0] = "Reaction"; 314 title[0] = "Reaction";
545 title[1] = "Reaction Rate [dm3/(mol*s)]"; 315 title[1] = "Reaction Rate [dm3/(mol*s)]";
>> 316 title[2] = "Interaction Range for chosen reaction model";
546 317
547 if (pReactionModel != nullptr) title[2] = << 318 G4cout<< setfill(' ')
548 "Interaction Range for chosen reaction << 319 << setw(maxlengthOutputReaction) << left << title[0]
549 << 320 << setw(maxlengthOutputReactionRate) << left << title[1]
550 G4cout << setfill(' ') << setw(maxlengthOu << 321 << setw(2) << left << title[2]
551 << setw(maxlengthOutputReactionRate) < << 322 << G4endl;
552 <<
553 if (pReactionModel != nullptr) G4cout << s <<
554 <<
555 G4cout << G4endl; <<
556 323
557 G4cout.fill('-'); 324 G4cout.fill('-');
558 if (pReactionModel != nullptr) G4cout.widt << 325 G4cout.width(maxlengthOutputReaction+2+maxlengthOutputReactionRate+2+(G4int)title[2].length());
559 maxlengthOutputReaction + 2 + maxlengt << 326 G4cout<<"-"<<G4endl;
560 + (G4int)title[2].length()); <<
561 else G4cout.width(maxlengthOutputReaction <<
562 G4cout << "-" << G4endl; <<
563 G4cout.fill(' '); 327 G4cout.fill(' ');
564 328
565 for (G4int i = 0; i < n; i++) << 329 for(G4int i = 0 ; i < n ; i ++)
566 { 330 {
567 G4cout << setw(maxlengthOutputReaction << 331 G4cout<< setw(maxlengthOutputReaction)<< left << outputReaction[i]
568 << setw(maxlengthOutputReaction << 332 << setw(maxlengthOutputReactionRate) << left << outputReactionRate[i]
569 << outputReactionRate[i]; << 333 << setw(2) << left <<outputRange[i]
570 << 334 <<G4endl;
571 if (pReactionModel != nullptr) G4cout <<
572 <<
573 G4cout << G4endl; <<
574 335
575 G4cout.fill('-'); 336 G4cout.fill('-');
576 if (pReactionModel != nullptr) G4cout. << 337 G4cout.width(maxlengthOutputReaction+2+maxlengthOutputReactionRate+2+(G4int)title[2].length());
577 maxlengthOutputReaction + 2 + maxl << 338 G4cout<<"-"<<G4endl;
578 + (G4int)title[2].length()); <<
579 else G4cout.width( <<
580 maxlengthOutputReaction + 2 + maxl <<
581 G4cout << "-" << G4endl; <<
582 G4cout.fill(' '); 339 G4cout.fill(' ');
583 } 340 }
584 341
585 delete[] title; << 342 delete [] outputReaction;
586 delete[] outputReaction; << 343 delete [] outputReactionRate;
587 delete[] outputReactionRate; << 344 delete [] outputRange;
588 delete[] outputRange; <<
589 } 345 }
590 << 346 //_____________________________________________________________________________________
591 //____________________________________________ <<
592 // Get/Set methods 347 // Get/Set methods
593 348
594 G4VDNAMolecularGeometry* G4DNAMolecularReactio << 349 const G4DNAMolecularReactionData*
>> 350 G4DNAMolecularReactionTable::GetReactionData(const G4Molecule* reactive1,
>> 351 const G4Molecule* reactive2) const
595 { 352 {
596 return fGeometry; << 353 if(fReactionData.empty())
597 } <<
598 <<
599 G4DNAMolecularReactionTable::Data* <<
600 G4DNAMolecularReactionTable::GetReactionData(R <<
601 R <<
602 { <<
603 if (fReactionData.empty()) <<
604 { 354 {
605 G4String errMsg = "No reaction table w 355 G4String errMsg = "No reaction table was implemented";
606 G4Exception("G4MolecularInteractionTab << 356 G4Exception("G4MolecularInteractionTable::CanInteractWith","",FatalErrorInArgument, errMsg);
607 FatalErrorInArgument, errM << 357 return 0;
608 } 358 }
609 359
610 auto it1 = fReactionData.find(pReactant1); << 360 ReactionDataMap::const_iterator it1 = fReactionData.find(reactive1);
611 361
612 if (it1 == fReactionData.end()) << 362 if(it1 == fReactionData.end())
613 { 363 {
614 G4String errMsg = << 364 G4cout<<"Nom : " << reactive1->GetName()<<G4endl;
615 "No reaction table was implemented << 365 G4String errMsg = "No reaction table was implemented for this molecule Definition : "
616 ->GetName(); << 366 + reactive1 -> GetName();
617 G4Exception("G4MolecularInteractionTab << 367 G4Exception("G4MolecularInteractionTable::CanReactWith","",FatalErrorInArgument, errMsg);
618 FatalErrorInArgument, errM <<
619 // Though the above is Fatal and will <<
620 return nullptr; <<
621 } 368 }
622 369
623 auto it2 = it1->second.find(pReactant2); << 370 std::map<const G4Molecule*,
>> 371 const G4DNAMolecularReactionData*,
>> 372 compMoleculeP>::const_iterator it2 = it1->second.find(reactive2);
624 373
625 if (it2 == it1->second.end()) << 374 if(it2 == it1->second.end())
626 { 375 {
627 G4cout << "Name : " << pReactant2->Get << 376 G4cout<<"Nom : " << reactive2->GetName()<<G4endl;
628 G4String errMsg = "No reaction table w << 377 G4String errMsg = "No reaction table was implemented for this molecule Definition : "
629 + pReactant2->GetName(); << 378 + reactive2 -> GetName();
630 G4Exception("G4MolecularInteractionTab << 379 G4Exception("G4MolecularInteractionTable::CanReactWith","",FatalErrorInArgument, errMsg);
631 } 380 }
632 381
633 return (it2->second); 382 return (it2->second);
634 } 383 }
635 384
636 const G4DNAMolecularReactionTable::ReactionDat << 385 const std::vector<const G4Molecule*>*
637 { << 386 G4DNAMolecularReactionTable::CanReactWith(const G4Molecule * aMolecule) const
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 { 387 {
658 if (fReactantsMV.empty()) << 388 if(fReactivesMV.empty())
659 { 389 {
660 G4String errMsg = "No reaction table w 390 G4String errMsg = "No reaction table was implemented";
661 G4Exception("G4MolecularInteractionTab << 391 G4Exception("G4MolecularInteractionTable::CanReactWith","",FatalErrorInArgument, errMsg);
662 FatalErrorInArgument, errM << 392 return 0;
663 return nullptr; <<
664 } 393 }
665 394
666 auto itReactivesMap = fReactantsMV.find(pM << 395 ReactivesMV::const_iterator itReactivesMap = fReactivesMV.find(aMolecule) ;
667 396
668 if (itReactivesMap == fReactantsMV.end()) << 397 if(itReactivesMap == fReactivesMV.end())
669 { 398 {
670 #ifdef G4VERBOSE << 399 G4cout<<"Nom : " << aMolecule->GetName()<<G4endl;
671 if (fVerbose) << 400 G4String errMsg = "No reaction table was implemented for this molecule Definition : "
672 { << 401 + aMolecule -> GetName();
673 G4String errMsg = "No reaction tab << 402 G4Exception("G4MolecularInteractionTable::CanReactWith","",FatalErrorInArgument, errMsg);
674 + pMolecule->GetName(); << 403 return 0;
675 // G4Exception("G4Molecul <<
676 G4cout << "--- G4MolecularInteract <<
677 G4cout << errMsg << G4endl; <<
678 } <<
679 #endif <<
680 return nullptr; <<
681 } 404 }
682 << 405 else
683 if (fVerbose) <<
684 { 406 {
685 G4cout << " G4MolecularInteractionTabl << 407 if(fVerbose)
686 G4cout << "You are checking reactants << 408 {
687 G4cout << " the number of reactants is << 409 G4cout<< " G4MolecularInteractionTable::CanReactWith :"<<G4endl;
>> 410 G4cout<<"You are checking reactants for : " << aMolecule->GetName()<<G4endl;
>> 411 G4cout<<" the number of reactants is : " << itReactivesMap->second.size()<<G4endl;
688 412
689 auto itProductsVector = itReactivesMap << 413 std::vector<const G4Molecule*>::const_iterator itProductsVector =
>> 414 itReactivesMap->second.begin();
690 415
691 for (; itProductsVector != itReactives << 416 for( ; itProductsVector != itReactivesMap->second.end(); itProductsVector++)
692 { << 417 {
693 G4cout << (*itProductsVector)->Get << 418 G4cout<<(*itProductsVector)->GetName()<<G4endl;
>> 419 }
694 } 420 }
>> 421 return &(itReactivesMap->second);
695 } 422 }
696 return &(itReactivesMap->second); << 423 return 0;
697 } 424 }
698 425
699 //____________________________________________ << 426 //_____________________________________________________________________________________
700 << 427 const std::map<const G4Molecule*, const G4DNAMolecularReactionData*, compMoleculeP>*
701 const G4DNAMolecularReactionTable::SpecificDat << 428 G4DNAMolecularReactionTable::GetReativesNData(const G4Molecule* molecule) const
702 G4DNAMolecularReactionTable::GetReativesNData( <<
703 { 429 {
704 if (fReactionData.empty()) << 430
>> 431 if(fReactionData.empty())
705 { 432 {
706 G4String errMsg = "No reaction table w 433 G4String errMsg = "No reaction table was implemented";
707 G4Exception("G4MolecularInteractionTab << 434 G4Exception("G4MolecularInteractionTable::CanInteractWith","",FatalErrorInArgument, errMsg);
708 FatalErrorInArgument, errM << 435 return 0;
709 } 436 }
710 437
711 auto itReactivesMap = fReactionData.find(m << 438 ReactionDataMap::const_iterator itReactivesMap = fReactionData.find(molecule) ;
712 439
713 if (itReactivesMap == fReactionData.end()) << 440 if(itReactivesMap == fReactionData.end())
714 { 441 {
715 return nullptr; << 442 G4cout<<"Nom : " << molecule->GetName()<<G4endl;
>> 443 G4String errMsg = "No reaction table was implemented for this molecule Definition : "
>> 444 + molecule -> GetName();
>> 445 G4Exception("G4MolecularInteractionTable::CanReactWith","",FatalErrorInArgument, errMsg);
716 } 446 }
717 << 447 else
718 if (fVerbose) <<
719 { 448 {
720 G4cout << " G4MolecularInteractionTabl << 449 if(fVerbose)
721 G4cout << "You are checking reactants <<
722 G4cout << " the number of reactants is <<
723 <<
724 auto itProductsVector = itReactivesMap <<
725 <<
726 for (; itProductsVector != itReactives <<
727 { 450 {
728 G4cout << itProductsVector->first- << 451 G4cout<< " G4MolecularInteractionTable::CanReactWith :"<<G4endl;
>> 452 G4cout<<"You are checking reactants for : " << molecule->GetName()<<G4endl;
>> 453 G4cout<<" the number of reactants is : " << itReactivesMap->second.size()<<G4endl;
>> 454
>> 455 std::map<const G4Molecule*,
>> 456 const G4DNAMolecularReactionData*,
>> 457 compMoleculeP>::const_iterator itProductsVector =
>> 458 itReactivesMap->second.begin();
>> 459
>> 460 for( ; itProductsVector != itReactivesMap->second.end(); itProductsVector++)
>> 461 {
>> 462 G4cout<<itProductsVector->first->GetName()<<G4endl;
>> 463 }
729 } 464 }
>> 465 return &(itReactivesMap->second);
730 } 466 }
731 return &(itReactivesMap->second); <<
732 } <<
733 467
734 //____________________________________________ << 468 return 0;
>> 469 }
735 470
736 const G4DNAMolecularReactionTable::DataList* << 471 const std::vector<const G4DNAMolecularReactionData*>*
737 G4DNAMolecularReactionTable::GetReactionData(c << 472 G4DNAMolecularReactionTable::GetReactionData(const G4Molecule* molecule) const
738 { 473 {
739 if (fReactionDataMV.empty()) << 474 if(fReactionDataMV.empty())
740 { 475 {
741 G4String errMsg = "No reaction table w 476 G4String errMsg = "No reaction table was implemented";
742 G4Exception("G4MolecularInteractionTab << 477 G4Exception("G4MolecularInteractionTable::CanInteractWith","",FatalErrorInArgument, errMsg);
743 FatalErrorInArgument, errM << 478 return 0 ;
744 } 479 }
745 auto it = fReactionDataMV.find(molecule); << 480 ReactionDataMV::const_iterator it = fReactionDataMV.find(molecule) ;
746 481
747 if (it == fReactionDataMV.end()) << 482 if(it == fReactionDataMV.end())
748 { 483 {
>> 484 G4cout<<"Nom : " << molecule->GetName()<<G4endl;
749 G4String errMsg = "No reaction table w 485 G4String errMsg = "No reaction table was implemented for this molecule Definition : "
750 + molecule->GetName(); << 486 + molecule -> GetName();
751 G4Exception("G4MolecularInteractionTab << 487 G4Exception("G4MolecularInteractionTable::GetReactionData","",FatalErrorInArgument, errMsg);
752 // Though the above is Fatal and will << 488 return 0; // coverity
753 return nullptr; <<
754 } 489 }
755 490
756 return &(it->second); 491 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 } 492 }
844 493