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