| 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 /* 26 /*
27 * G4ReactionTableMessenger.cc 27 * G4ReactionTableMessenger.cc
28 * 28 *
29 * Created on: Sep 14, 2015 29 * Created on: Sep 14, 2015
30 * Author: mkaramit 30 * Author: mkaramit
31 */ 31 */
32 32
33 #include <G4ReactionTableMessenger.hh> 33 #include <G4ReactionTableMessenger.hh>
34 #include <G4DNAMolecularReactionTable.hh> 34 #include <G4DNAMolecularReactionTable.hh>
35 #include <G4UIcmdWithAString.hh> 35 #include <G4UIcmdWithAString.hh>
36 #include <G4UIcmdWithADoubleAndUnit.hh> 36 #include <G4UIcmdWithADoubleAndUnit.hh>
37 #include <G4UnitsTable.hh> 37 #include <G4UnitsTable.hh>
38 #include <G4SystemOfUnits.hh> 38 #include <G4SystemOfUnits.hh>
39 #include <G4UIcmdWithoutParameter.hh> 39 #include <G4UIcmdWithoutParameter.hh>
40 #include <G4UIcmdWithABool.hh> <<
41 40
42 //-------------------------------------------- 41 //------------------------------------------------------------------------------
43 42
44 G4ReactionTableMessenger::G4ReactionTableMesse << 43 G4ReactionTableMessenger::G4ReactionTableMessenger(G4DNAMolecularReactionTable* table) :
45 : << 44 G4UImessenger()
46 fpTable(table) <<
47 , fpActivateReactionUI(new G4UIcmdWithoutPar <<
48 { 45 {
>> 46 fpTable = table;
>> 47
49 fpNewDiffContReaction = new G4UIcmdWithAStri 48 fpNewDiffContReaction = new G4UIcmdWithAString("/chem/reaction/new", this);
50 fpAddReaction = new G4UIcmdWithAString("/che 49 fpAddReaction = new G4UIcmdWithAString("/chem/reaction/add", this);
>> 50 // fpNewPartDiffContReactionByRadius =
>> 51 // new G4UIcmdWithAString("/chem/reaction/newPartDiffByRadius", this);
>> 52 //
>> 53 // fpNewPartDiffContReactionByReactionRate =
>> 54 // new G4UIcmdWithAString("/chem/reaction/newPartDiffByRate", this);
>> 55
51 fpPrintTable = new G4UIcmdWithoutParameter(" 56 fpPrintTable = new G4UIcmdWithoutParameter("/chem/reaction/print", this);
52 } 57 }
53 58
54 //-------------------------------------------- 59 //------------------------------------------------------------------------------
55 60
56 G4ReactionTableMessenger::~G4ReactionTableMess 61 G4ReactionTableMessenger::~G4ReactionTableMessenger()
57 { 62 {
58 delete fpNewDiffContReaction; << 63 if(fpNewDiffContReaction) delete fpNewDiffContReaction;
59 delete fpAddReaction; << 64 if(fpAddReaction) delete fpAddReaction;
60 delete fpPrintTable; << 65 if(fpPrintTable) delete fpPrintTable;
61 } 66 }
62 67
63 //-------------------------------------------- 68 //------------------------------------------------------------------------------
>> 69
64 void G4ReactionTableMessenger::SetNewValue(G4U 70 void G4ReactionTableMessenger::SetNewValue(G4UIcommand* command,
65 G4S 71 G4String newValue)
66 { 72 {
67 if(command == fpActivateReactionUI.get()) <<
68 { <<
69 //assert(false); <<
70 fpTable->Reset();//release reaction data <<
71 } <<
72 73
73 if(command == fpNewDiffContReaction) 74 if(command == fpNewDiffContReaction)
74 { 75 {
75 std::istringstream iss(newValue); 76 std::istringstream iss(newValue);
76 77
77 G4String species1; 78 G4String species1;
78 iss >> species1; 79 iss >> species1;
79 80
80 G4String species2; 81 G4String species2;
81 iss >> species2; 82 iss >> species2;
82 83
83 double reactionRate; 84 double reactionRate;
84 iss >> reactionRate; 85 iss >> reactionRate;
85 86
86 // G4String reactionRateUnit; 87 // G4String reactionRateUnit;
87 // iss >> reactionRateUnit; 88 // iss >> reactionRateUnit;
88 89
89 double dimensionedReactionRate = reactionR 90 double dimensionedReactionRate = reactionRate * (1e-3 * m3 / (mole * s));
90 // G4UIcmdWithADoubleAndUnit::ConvertTo 91 // G4UIcmdWithADoubleAndUnit::ConvertToDimensionedDouble((reactionRate
91 // + G4String(" ") + reactionRateUn 92 // + G4String(" ") + reactionRateUnit).c_str());
92 93
93 auto reactionData = << 94 G4DNAMolecularReactionData* reactionData =
94 new G4DNAMolecularReactionData(dimensi 95 new G4DNAMolecularReactionData(dimensionedReactionRate,
95 species 96 species1,
96 species 97 species2);
97 98
98 // G4String productionRate; 99 // G4String productionRate;
99 // iss >> productionRate; 100 // iss >> productionRate;
100 // 101 //
101 // if(productionRate != "" && productionRat 102 // if(productionRate != "" && productionRate != "X")
102 // { 103 // {
103 // double prodRateReal = G4UIcommand::Con 104 // double prodRateReal = G4UIcommand::ConvertToDouble(productionRate);
104 // 105 //
105 // if(prodRateReal == 0 || isnan(prodRate 106 // if(prodRateReal == 0 || isnan(prodRateReal))
106 // { 107 // {
107 // G4Exception("G4ReactionTableMessenge 108 // G4Exception("G4ReactionTableMessenger",
108 // "WRONG_PRODUCTION_RATE", 109 // "WRONG_PRODUCTION_RATE",
109 // FatalException, ""); 110 // FatalException, "");
110 // } 111 // }
111 // 112 //
112 // double dimensionedProductionRate = pro 113 // double dimensionedProductionRate = prodRateReal
113 // * (1e-3 * m3 / (mole * s)); 114 // * (1e-3 * m3 / (mole * s));
114 // reactionData->SetProductionRate(dimens 115 // reactionData->SetProductionRate(dimensionedProductionRate);
115 // } 116 // }
116 117
117 while(!iss.eof()) << 118 while(iss.eof() == false)
118 { 119 {
119 G4String product; 120 G4String product;
120 iss >> product; 121 iss >> product;
121 122
122 if(!product.empty()) << 123 if(product != "")
123 { 124 {
124 reactionData->AddProduct(product); 125 reactionData->AddProduct(product);
125 } 126 }
126 else 127 else
127 { 128 {
128 break; 129 break;
129 } 130 }
130 }; 131 };
131 132
132 fpTable->SetReaction(reactionData); 133 fpTable->SetReaction(reactionData);
133 } 134 }
134 // else if(command == fpNewPartDiffContReacti 135 // else if(command == fpNewPartDiffContReactionByRadius)
135 // { 136 // {
136 // std::istringstream iss(newValue); 137 // std::istringstream iss(newValue);
137 // 138 //
138 // G4String species1; 139 // G4String species1;
139 // iss >> species1; 140 // iss >> species1;
140 // 141 //
141 // G4String species2; 142 // G4String species2;
142 // iss >> species2; 143 // iss >> species2;
143 // 144 //
144 // double reactionRate; 145 // double reactionRate;
145 // iss >> reactionRate; 146 // iss >> reactionRate;
146 // 147 //
147 //// G4String reactionRateUnit; 148 //// G4String reactionRateUnit;
148 //// iss >> reactionRateUnit; 149 //// iss >> reactionRateUnit;
149 // 150 //
150 //// G4String reactionRateUnit; 151 //// G4String reactionRateUnit;
151 //// iss >> reactionRateUnit; 152 //// iss >> reactionRateUnit;
152 // 153 //
153 // double reactionRadius; 154 // double reactionRadius;
154 // iss >> reactionRadius; 155 // iss >> reactionRadius;
155 // 156 //
156 //// G4String reactionRadiusUnit; 157 //// G4String reactionRadiusUnit;
157 //// iss >> reactionRadiusUnit; 158 //// iss >> reactionRadiusUnit;
158 // 159 //
159 // double dimensionedReactionRate = reactio 160 // double dimensionedReactionRate = reactionRate * (1e-3 * m3 / (mole * s));
160 // 161 //
161 //// double dimensionedReactionRate = 162 //// double dimensionedReactionRate =
162 //// G4UIcmdWithADoubleAndUnit::Convert 163 //// G4UIcmdWithADoubleAndUnit::ConvertToDimensionedDouble((reactionRate
163 //// + " " + reactionRateUnit).c_st 164 //// + " " + reactionRateUnit).c_str());
164 // 165 //
165 // double dimensionedReactionRadius = react 166 // double dimensionedReactionRadius = reactionRadius * nm;
166 //// G4UIcmdWithADoubleAndUnit::Convert 167 //// G4UIcmdWithADoubleAndUnit::ConvertToDimensionedDouble((reactionRadius
167 //// + " " + reactionRadiusUnit).c_ 168 //// + " " + reactionRadiusUnit).c_str());
168 // 169 //
169 // G4DNAMolecularReactionData* reactionData 170 // G4DNAMolecularReactionData* reactionData =
170 // new G4DNAMolecularReactionData(dimen 171 // new G4DNAMolecularReactionData(dimensionedReactionRate,
171 // speci 172 // species1,
172 // speci 173 // species2);
173 // reactionData->SetPartiallyDiffusionContr 174 // reactionData->SetPartiallyDiffusionControlledReaction(dimensionedReactionRate,
174 // 175 // dimensionedReactionRadius);
175 // 176 //
176 // while(iss.eof() == false) 177 // while(iss.eof() == false)
177 // { 178 // {
178 // G4String product; 179 // G4String product;
179 // iss >> product; 180 // iss >> product;
180 // 181 //
181 // if(product != "") 182 // if(product != "")
182 // { 183 // {
183 // reactionData->AddProduct(product); 184 // reactionData->AddProduct(product);
184 // } 185 // }
185 // else 186 // else
186 // { 187 // {
187 // break; 188 // break;
188 // } 189 // }
189 // } 190 // }
190 // 191 //
191 // fpTable->SetReaction(reactionData); 192 // fpTable->SetReaction(reactionData);
192 // } 193 // }
193 // else if(command == fpNewPartDiffContReacti 194 // else if(command == fpNewPartDiffContReactionByReactionRate)
194 // { 195 // {
195 // std::istringstream iss(newValue); 196 // std::istringstream iss(newValue);
196 // 197 //
197 // G4String species1; 198 // G4String species1;
198 // iss >> species1; 199 // iss >> species1;
199 // 200 //
200 // G4String species2; 201 // G4String species2;
201 // iss >> species2; 202 // iss >> species2;
202 // 203 //
203 // double reactionRate; 204 // double reactionRate;
204 // iss >> reactionRate; 205 // iss >> reactionRate;
205 // 206 //
206 // // G4String reactionRateUnit; 207 // // G4String reactionRateUnit;
207 // // iss >> reactionRateUnit; 208 // // iss >> reactionRateUnit;
208 // 209 //
209 // // G4String reactionRateUnit; 210 // // G4String reactionRateUnit;
210 // // iss >> reactionRateUnit; 211 // // iss >> reactionRateUnit;
211 // 212 //
212 // double activationRate; 213 // double activationRate;
213 // iss >> activationRate; 214 // iss >> activationRate;
214 // 215 //
215 // // G4String reactionRadiusUnit; 216 // // G4String reactionRadiusUnit;
216 // // iss >> reactionRadiusUnit; 217 // // iss >> reactionRadiusUnit;
217 // 218 //
218 // double dimensionedReactionRate = reactio 219 // double dimensionedReactionRate = reactionRate * (1e-3 * m3 / (mole * s));
219 // 220 //
220 // double dimensionedActivationRate = activ 221 // double dimensionedActivationRate = activationRate
221 // * (1e-3 * m3 / (mole * s)); 222 // * (1e-3 * m3 / (mole * s));
222 // 223 //
223 // G4DNAMolecularReactionData* reactionData 224 // G4DNAMolecularReactionData* reactionData =
224 // new G4DNAMolecularReactionData(dimen 225 // new G4DNAMolecularReactionData(dimensionedReactionRate,
225 // speci 226 // species1,
226 // speci 227 // species2);
227 // reactionData-> 228 // reactionData->
228 // SetPartiallyDiffusionControlledReactio 229 // SetPartiallyDiffusionControlledReactionByActivation(dimensionedReactionRate,
229 // 230 // dimensionedActivationRate);
230 // 231 //
231 // while(iss.eof() == false) 232 // while(iss.eof() == false)
232 // { 233 // {
233 // G4String product; 234 // G4String product;
234 // iss >> product; 235 // iss >> product;
235 // 236 //
236 // if(product != "") 237 // if(product != "")
237 // { 238 // {
238 // reactionData->AddProduct(product); 239 // reactionData->AddProduct(product);
239 // } 240 // }
240 // else 241 // else
241 // { 242 // {
242 // break; 243 // break;
243 // } 244 // }
244 // } 245 // }
245 // 246 //
246 // fpTable->SetReaction(reactionData); 247 // fpTable->SetReaction(reactionData);
247 // } 248 // }
248 else if(command == fpPrintTable) 249 else if(command == fpPrintTable)
249 { 250 {
250 fpTable->PrintTable(); 251 fpTable->PrintTable();
251 } 252 }
252 else if(command == fpAddReaction) 253 else if(command == fpAddReaction)
253 { 254 {
254 std::istringstream iss(newValue); 255 std::istringstream iss(newValue);
255 256
256 //---------------------------------------- 257 //--------------------------------------------------------------------------
257 // Reactants definition 258 // Reactants definition
258 259
259 G4String species1; 260 G4String species1;
260 iss >> species1; 261 iss >> species1;
261 262
262 G4String marker; 263 G4String marker;
263 iss >> marker; // peut etre +, ->, | 264 iss >> marker; // peut etre +, ->, |
264 265
265 G4String species2; 266 G4String species2;
266 267
267 if(marker == "+") 268 if(marker == "+")
268 { 269 {
269 iss >> species2; 270 iss >> species2;
270 iss >> marker; // peut etre ->, | 271 iss >> marker; // peut etre ->, |
271 } 272 }
272 273
273 //---------------------------------------- 274 //--------------------------------------------------------------------------
274 275
275 auto reactionData = << 276 G4DNAMolecularReactionData* reactionData =
276 new G4DNAMolecularReactionData 277 new G4DNAMolecularReactionData(0,
277 278 species1,
278 279 species2);
279 //fpTable->SetReaction(reactionData); 280 //fpTable->SetReaction(reactionData);
280 281
281 //---------------------------------------- 282 //--------------------------------------------------------------------------
282 // Add products 283 // Add products
283 if(marker == "->") 284 if(marker == "->")
284 { 285 {
285 iss >> marker; // doit etre = species na 286 iss >> marker; // doit etre = species name
286 287
287 while(marker!="|" 288 while(marker!="|"
288 //&& marker!="" 289 //&& marker!=""
289 && !iss.eof() << 290 && iss.eof() == false
290 ) 291 )
291 { 292 {
292 //G4cout << marker << G4endl; << 293 G4cout << marker << G4endl;
293 if(marker == "+") 294 if(marker == "+")
294 { 295 {
295 iss >> marker; // doit etre species 296 iss >> marker; // doit etre species name
296 continue; 297 continue;
297 } 298 }
298 if(marker != "H2O") << 299 reactionData->AddProduct(marker);
299 { <<
300 reactionData->AddProduct(marker); <<
301 } <<
302 <<
303 iss >> marker; // peut etre species na 300 iss >> marker; // peut etre species name, +, |
304 }; 301 };
305 } 302 }
306 303
307 // G4cout << "out of 1st loop" << G4endl; 304 // G4cout << "out of 1st loop" << G4endl;
308 305
309 //---------------------------------------- 306 //--------------------------------------------------------------------------
310 // Add reaction rate method 307 // Add reaction rate method
311 G4String rateconst_method; 308 G4String rateconst_method;
312 iss >> rateconst_method; 309 iss >> rateconst_method;
313 if(rateconst_method == "Fix") 310 if(rateconst_method == "Fix")
314 { 311 {
315 iss >> marker; // must be | 312 iss >> marker; // must be |
316 double reactionRate; 313 double reactionRate;
317 iss >> reactionRate; 314 iss >> reactionRate;
318 315
319 double dimensionedReactionRate = reactio 316 double dimensionedReactionRate = reactionRate * (1e-3 * m3 / (mole * s));
320 reactionData->SetObservedReactionRateCon 317 reactionData->SetObservedReactionRateConstant(dimensionedReactionRate);
321 reactionData->ComputeEffectiveRadius(); <<
322 G4String markerType; <<
323 iss >> markerType; // must be | <<
324 if(markerType == "|") <<
325 { <<
326 G4int reactionType; <<
327 iss >> reactionType; <<
328 if(reactionType == 1) <<
329 { <<
330 reactionData->SetReactionType(reacti <<
331 } <<
332 } <<
333 <<
334 318
335 // G4String productionRate; 319 // G4String productionRate;
336 // iss >> productionRate; 320 // iss >> productionRate;
337 // 321 //
338 // if(productionRate != "" && productionR 322 // if(productionRate != "" && productionRate != "X")
339 // { 323 // {
340 // double prodRateReal = G4UIcommand::C 324 // double prodRateReal = G4UIcommand::ConvertToDouble(productionRate);
341 // 325 //
342 // if(prodRateReal == 0 || isnan(prodRa 326 // if(prodRateReal == 0 || isnan(prodRateReal))
343 // { 327 // {
344 // G4Exception("G4ReactionTableMessen 328 // G4Exception("G4ReactionTableMessenger",
345 // "WRONG_PRODUCTION_RATE 329 // "WRONG_PRODUCTION_RATE",
346 // FatalException, 330 // FatalException,
347 // ""); 331 // "");
348 // } 332 // }
349 // 333 //
350 // double dimensionedProductionRate = p 334 // double dimensionedProductionRate = prodRateReal
351 // * (1e-3 * m3 / (mole * s)); 335 // * (1e-3 * m3 / (mole * s));
352 // reactionData->SetProductionRate(dime 336 // reactionData->SetProductionRate(dimensionedProductionRate);
353 // } 337 // }
354 } 338 }
355 else if(rateconst_method == "Arr") 339 else if(rateconst_method == "Arr")
356 { 340 {
357 iss >> marker; // must be | 341 iss >> marker; // must be |
358 double A0 = 0; 342 double A0 = 0;
359 double E_R = 0; 343 double E_R = 0;
360 344
361 iss >> A0; 345 iss >> A0;
362 iss >> E_R; 346 iss >> E_R;
363 reactionData->SetArrehniusParameterizati 347 reactionData->SetArrehniusParameterization(A0, E_R);
364 } 348 }
365 else if(rateconst_method == "Pol") 349 else if(rateconst_method == "Pol")
366 { 350 {
367 iss >> marker; // must be | 351 iss >> marker; // must be |
368 std::vector<double> P = {0, 0, 0, 0, 0}; 352 std::vector<double> P = {0, 0, 0, 0, 0};
369 353
370 size_t i = 0; 354 size_t i = 0;
371 while(i < 4) // could be changed to 5 on 355 while(i < 4) // could be changed to 5 only if marker is used as delimiter
372 { 356 {
373 double tmp; 357 double tmp;
374 iss >> tmp; 358 iss >> tmp;
375 P[i] = tmp; 359 P[i] = tmp;
376 // G4cout << newValue << G4endl; 360 // G4cout << newValue << G4endl;
377 // G4cout << tmp << G4endl; 361 // G4cout << tmp << G4endl;
378 // G4cout << P[i] << G4endl; 362 // G4cout << P[i] << G4endl;
379 ++i; 363 ++i;
380 }; 364 };
381 reactionData->SetPolynomialParameterizat 365 reactionData->SetPolynomialParameterization(P);
382 } 366 }
383 else if(rateconst_method == "Scale") 367 else if(rateconst_method == "Scale")
384 { 368 {
385 iss >> marker; // must be | 369 iss >> marker; // must be |
386 double temp_K; 370 double temp_K;
387 iss >> temp_K; 371 iss >> temp_K;
388 double reactionRateCste; 372 double reactionRateCste;
389 iss >> reactionRateCste; 373 iss >> reactionRateCste;
390 double dimensionedReactionRate = reactio 374 double dimensionedReactionRate = reactionRateCste * (1e-3 * m3 / (mole * s));
391 reactionData->SetObservedReactionRateCon 375 reactionData->SetObservedReactionRateConstant(dimensionedReactionRate);
392 reactionData->SetScaledParameterization( 376 reactionData->SetScaledParameterization(temp_K, dimensionedReactionRate);
393 } 377 }
394 378
395 // if(iss.eof() == false) 379 // if(iss.eof() == false)
396 // { 380 // {
397 // iss >> marker; 381 // iss >> marker;
398 // 382 //
399 // if(marker == "|") 383 // if(marker == "|")
400 // { 384 // {
401 // G4String productionRate ; 385 // G4String productionRate ;
402 // iss >> productionRate; 386 // iss >> productionRate;
403 // 387 //
404 //// G4cout << productionRate << G4endl 388 //// G4cout << productionRate << G4endl;
405 // 389 //
406 // double dimProductionRate = G4UIcomma 390 // double dimProductionRate = G4UIcommand::ConvertToDouble(productionRate)* (1e-3 * m3 / (mole * s));
407 // 391 //
408 // G4cout << " DIM PROD RATE = " << rea 392 // G4cout << " DIM PROD RATE = " << reactionData->GetReactant1()->GetName()
409 // << " + " << reactionData->GetRe 393 // << " + " << reactionData->GetReactant2()->GetName() << " = " << dimProductionRate << G4endl;
410 // 394 //
411 // reactionData->SetProductionRate(dimP 395 // reactionData->SetProductionRate(dimProductionRate);
412 // } 396 // }
413 // } 397 // }
414 fpTable->SetReaction(reactionData); 398 fpTable->SetReaction(reactionData);
415 // G4cout << "Reaction " << species1 << " + 399 // G4cout << "Reaction " << species1 << " + " << species2 << " added" << G4endl;
416 } 400 }
417 } 401 }
418 402