| 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 #include <memory> 27 #include <memory>
28 #include "G4DNAEventScheduler.hh" 28 #include "G4DNAEventScheduler.hh"
29 #include "G4DNAGillespieDirectMethod.hh" 29 #include "G4DNAGillespieDirectMethod.hh"
30 #include "G4SystemOfUnits.hh" 30 #include "G4SystemOfUnits.hh"
31 #include "G4UnitsTable.hh" 31 #include "G4UnitsTable.hh"
32 #include "G4DNAUpdateSystemModel.hh" 32 #include "G4DNAUpdateSystemModel.hh"
33 #include "G4DNAMolecularReactionTable.hh" 33 #include "G4DNAMolecularReactionTable.hh"
34 #include "G4Timer.hh" 34 #include "G4Timer.hh"
35 #include "G4Scheduler.hh" 35 #include "G4Scheduler.hh"
36 #include "G4UserMeshAction.hh" 36 #include "G4UserMeshAction.hh"
37 #include "G4MoleculeCounter.hh" 37 #include "G4MoleculeCounter.hh"
38 #include "G4DNAScavengerMaterial.hh" 38 #include "G4DNAScavengerMaterial.hh"
39 #include "G4Molecule.hh" <<
40 39
41 G4DNAEventScheduler::G4DNAEventScheduler() << 40 G4DNAEventScheduler::G4DNAEventScheduler(const G4DNABoundingBox& boundingBox,
42 : fpGillespieReaction(new G4DNAGillespieDire << 41 G4int pixel)
>> 42 : IEventScheduler()
>> 43 , fVerbose(0)
>> 44 , fInitialized(false)
>> 45 , fStartTime(1 * ps)
>> 46 , fEndTime(10000 * s)
>> 47 , fStepNumber(0)
>> 48 , fMaxStep(INT_MAX)
>> 49 , fRunning(true)
>> 50 , fTimeStep(DBL_MAX)
>> 51 , fGlobalTime(fStartTime)
>> 52 , fJumpingNumber(0)
>> 53 , fReactionNumber(0)
>> 54 , fPixel(pixel)
>> 55 , fIsChangeMesh(false)
>> 56 , fSetChangeMesh(true)
>> 57 , fStepNumberInMesh(0)
>> 58 , fInitialPixels(fPixel)
>> 59 , fpMesh(new G4DNAMesh(boundingBox, fPixel))
>> 60 , fpGillespieReaction(new G4DNAGillespieDirectMethod())
43 , fpEventSet(new G4DNAEventSet()) 61 , fpEventSet(new G4DNAEventSet())
44 , fpUpdateSystem(new G4DNAUpdateSystemModel( 62 , fpUpdateSystem(new G4DNAUpdateSystemModel())
45 {} << 63 , fpUserMeshAction(nullptr)
>> 64 {
>> 65 if(!CheckingReactionRadius(fpMesh->GetResolution()))
>> 66 {
>> 67 G4String WarMessage = "resolution is not good : " +
>> 68 std::to_string(fpMesh->GetResolution() / nm);
>> 69 G4Exception("G4DNAEventScheduler::InitializeInMesh()", "WrongResolution",
>> 70 JustWarning, WarMessage);
>> 71 }
>> 72 }
46 73
47 void G4DNAEventScheduler::ClearAndReChargeCoun 74 void G4DNAEventScheduler::ClearAndReChargeCounter()
48 { 75 {
49 fCounterMap.clear(); 76 fCounterMap.clear();
50 if(fTimeToRecord.empty()) 77 if(fTimeToRecord.empty())
51 { 78 {
52 G4String WarMessage = "fTimeToRecord is em << 79 G4cout << "fTimeToRecord is empty " << G4endl;
53 G4Exception("G4DNAEventScheduler::ClearAnd <<
54 "TimeToRecord is empty", JustW <<
55 } 80 }
56 fLastRecoredTime = fTimeToRecord.begin(); 81 fLastRecoredTime = fTimeToRecord.begin();
57 82
58 if(G4VMoleculeCounter::Instance()->InUse()) 83 if(G4VMoleculeCounter::Instance()->InUse()) // copy from MoleculeCounter
59 { 84 {
60 G4MoleculeCounter::RecordedMolecules speci 85 G4MoleculeCounter::RecordedMolecules species;
61 species = G4MoleculeCounter::Instance()->G 86 species = G4MoleculeCounter::Instance()->GetRecordedMolecules();
62 if(species.get() == nullptr) 87 if(species.get() == nullptr)
63 { 88 {
64 return; 89 return;
65 } 90 }
66 if(species->empty()) << 91 else if(species->empty())
67 { 92 {
68 G4MoleculeCounter::Instance()->ResetCoun 93 G4MoleculeCounter::Instance()->ResetCounter();
69 return; 94 return;
70 } 95 }
71 for(auto time_mol : fTimeToRecord) 96 for(auto time_mol : fTimeToRecord)
72 { 97 {
73 if(time_mol > fStartTime) 98 if(time_mol > fStartTime)
74 { 99 {
75 continue; 100 continue;
76 } 101 }
77 102
78 for(auto molecule : *species) 103 for(auto molecule : *species)
79 { 104 {
80 G4int n_mol = G4MoleculeCounter::Insta 105 G4int n_mol = G4MoleculeCounter::Instance()->GetNMoleculesAtTime(
81 molecule, time_mol); 106 molecule, time_mol);
82 107
83 if(n_mol < 0) 108 if(n_mol < 0)
84 { 109 {
85 G4cerr << "G4DNAEventScheduler::Clea 110 G4cerr << "G4DNAEventScheduler::ClearAndReChargeCounter() ::N "
86 "molecules not valid < 0 " 111 "molecules not valid < 0 "
87 << G4endl; 112 << G4endl;
88 G4Exception("", "N<0", FatalExceptio 113 G4Exception("", "N<0", FatalException, "");
89 } 114 }
90 fCounterMap[time_mol][molecule] = n_mo 115 fCounterMap[time_mol][molecule] = n_mol;
91 } 116 }
>> 117
92 fLastRecoredTime++; 118 fLastRecoredTime++;
93 } 119 }
94 G4MoleculeCounter::Instance()->ResetCounte 120 G4MoleculeCounter::Instance()->ResetCounter(); // reset
95 G4MoleculeCounter::Instance()->Use(false); 121 G4MoleculeCounter::Instance()->Use(false); // no more used
96 } 122 }
97 else <<
98 { <<
99 G4ExceptionDescription exceptionDescriptio <<
100 exceptionDescription << "G4VMoleculeCounte <<
101 G4Exception("G4DNAEventScheduler::ClearAnd <<
102 "G4DNAEventScheduler010", Just <<
103 } <<
104 } 123 }
105 124
106 [[maybe_unused]] void G4DNAEventScheduler::Add 125 [[maybe_unused]] void G4DNAEventScheduler::AddTimeToRecord(const G4double& time)
107 { 126 {
108 if(fTimeToRecord.find(time) == fTimeToRecord 127 if(fTimeToRecord.find(time) == fTimeToRecord.end())
109 { 128 {
110 fTimeToRecord.insert(time); 129 fTimeToRecord.insert(time);
111 } 130 }
112 } 131 }
113 132
114 G4DNAEventScheduler::~G4DNAEventScheduler() = 133 G4DNAEventScheduler::~G4DNAEventScheduler() = default;
115 134
116 void G4DNAEventScheduler::Voxelizing() 135 void G4DNAEventScheduler::Voxelizing()
117 { 136 {
118 auto pMainList = G4ITTrackHolder::Instance() 137 auto pMainList = G4ITTrackHolder::Instance()->GetMainList();
119 std::map<G4VDNAMesh::Index, MapList> TrackKe << 138 std::map<G4DNAMesh::Key, MapList> TrackKeyMap;
120 for(auto track : *pMainList) 139 for(auto track : *pMainList)
121 { 140 {
122 auto molType = GetMolecule(track)->GetMole 141 auto molType = GetMolecule(track)->GetMolecularConfiguration();
123 142
124 auto pScavengerMaterial = dynamic_cast<G4D 143 auto pScavengerMaterial = dynamic_cast<G4DNAScavengerMaterial*>(
125 G4Scheduler::Instance()->GetScavengerMat 144 G4Scheduler::Instance()->GetScavengerMaterial());
126 if(pScavengerMaterial != nullptr && 145 if(pScavengerMaterial != nullptr &&
127 pScavengerMaterial->find(molType)) // 146 pScavengerMaterial->find(molType)) // avoid voxelize the scavenger
128 { 147 {
129 continue; 148 continue;
130 } 149 }
131 150
132 auto key = fpMesh->GetIndex(track->GetPosi << 151 auto key = fpMesh->GetKey(track->GetPosition());
133 if(TrackKeyMap.find(key) != TrackKeyMap.en 152 if(TrackKeyMap.find(key) != TrackKeyMap.end())
134 { 153 {
135 std::map<MolType, size_t>& TrackTypeMap 154 std::map<MolType, size_t>& TrackTypeMap = TrackKeyMap[key];
136 if(TrackTypeMap.find(molType) != TrackTy 155 if(TrackTypeMap.find(molType) != TrackTypeMap.end())
137 { 156 {
138 TrackTypeMap[molType]++; 157 TrackTypeMap[molType]++;
139 } 158 }
140 else 159 else
141 { 160 {
142 TrackTypeMap[molType] = 1; 161 TrackTypeMap[molType] = 1;
143 } 162 }
144 } 163 }
145 else 164 else
146 { 165 {
147 TrackKeyMap[key][molType] = 1; 166 TrackKeyMap[key][molType] = 1;
148 } 167 }
149 } 168 }
150 169
151 for(auto& it : TrackKeyMap) 170 for(auto& it : TrackKeyMap)
152 { 171 {
153 fpMesh->InitializeVoxel(it.first, std::mov << 172 fpMesh->SetVoxelMapList(it.first, std::move(it.second));
154 } 173 }
155 } 174 }
156 175
157 void G4DNAEventScheduler::ReVoxelizing(G4int p 176 void G4DNAEventScheduler::ReVoxelizing(G4int pixel)
158 { 177 {
159 fPixel = pixel; 178 fPixel = pixel;
160 auto newMesh = new G4DNAMesh(fpMesh->GetBoun 179 auto newMesh = new G4DNAMesh(fpMesh->GetBoundingBox(), fPixel);
161 180
162 auto begin = fpMesh->begin(); 181 auto begin = fpMesh->begin();
163 auto end = fpMesh->end(); 182 auto end = fpMesh->end();
164 std::map<G4VDNAMesh::Index, MapList> TrackKe << 183 std::map<G4DNAMesh::Key, MapList> TrackKeyMap;
165 for(; begin != end; begin++) 184 for(; begin != end; begin++)
166 { 185 {
167 auto index = std::get<0>(*begin); << 186 auto index = fpMesh->GetIndex(begin->first);
168 auto newIndex = fpMesh->ConvertIndex(index << 187 auto newKey = newMesh->GetKey(fpMesh->GetIndex(index, fPixel));
169 if(TrackKeyMap.find(newIndex) == TrackKeyM << 188 auto node = begin->second;
>> 189 // if (node == nullptr) continue;
>> 190 if(TrackKeyMap.find(newKey) == TrackKeyMap.end())
170 { 191 {
171 TrackKeyMap[newIndex] = std::get<2>(*beg << 192 TrackKeyMap[newKey] = node->GetMapList();
172 } 193 }
173 else 194 else
174 { 195 {
175 for(const auto& it : std::get<2>(*begin) << 196 for(const auto& it : node->GetMapList())
176 { 197 {
177 TrackKeyMap[newIndex][it.first] += it. << 198 TrackKeyMap[newKey][it.first] += it.second;
178 } 199 }
179 if(fVerbose > 1) 200 if(fVerbose > 1)
180 { 201 {
181 G4cout << " ReVoxelizing:: Old index : << 202 G4cout << "key : " << begin->first << " index : " << index
182 << " new index : " << fpMesh->C << 203 << " new index : " << fpMesh->GetIndex(index, fPixel)
183 << " number: " << std::get<2>(* << 204 << " new key : " << newKey
>> 205 << " number: " << node->GetMapList().begin()->second << G4endl;
184 } 206 }
185 } 207 }
186 } 208 }
187 fpMesh.reset(newMesh); 209 fpMesh.reset(newMesh);
188 210
189 for(auto& it : TrackKeyMap) 211 for(auto& it : TrackKeyMap)
190 { 212 {
191 fpMesh->InitializeVoxel(it.first, std::mov << 213 fpMesh->SetVoxelMapList(it.first, std::move(it.second));
192 } 214 }
193 } 215 }
194 void G4DNAEventScheduler::Reset() 216 void G4DNAEventScheduler::Reset()
195 { 217 {
196 // find another solultion 218 // find another solultion
197 fGlobalTime = fEndTime; 219 fGlobalTime = fEndTime;
198 220
199 // 221 //
200 LastRegisterForCounter();//Last register for << 222 // RecordTime();//Last register for counter
201 223
202 if(fVerbose > 0) 224 if(fVerbose > 0)
203 { 225 {
204 G4cout << "End Processing and reset Gird, 226 G4cout << "End Processing and reset Gird, ScavengerTable, EventSet for new "
205 "simulation!!!!" 227 "simulation!!!!"
206 << G4endl; 228 << G4endl;
207 } 229 }
208 fInitialized = false; 230 fInitialized = false;
209 fTimeStep = 0; 231 fTimeStep = 0;
210 fStepNumber = 0; 232 fStepNumber = 0;
211 fGlobalTime = fStartTime; 233 fGlobalTime = fStartTime;
212 fRunning = true; 234 fRunning = true;
213 fReactionNumber = 0; 235 fReactionNumber = 0;
214 fJumpingNumber = 0; 236 fJumpingNumber = 0;
215 237
216 fpEventSet->RemoveEventSet(); 238 fpEventSet->RemoveEventSet();
217 fpMesh->Reset(); 239 fpMesh->Reset();
218 fpGillespieReaction->ResetEquilibrium(); <<
219 } 240 }
220 241
221 void G4DNAEventScheduler::Initialize(const G4D << 242 void G4DNAEventScheduler::Initialize()
222 G4int pix <<
223 { 243 {
224 if(!fInitialized) 244 if(!fInitialized)
225 { 245 {
226 fPixel = pixel; << 246 fPixel = fInitialPixels;
227 fpMesh = std::make_unique<G4DNAMesh>(bound << 247 fpMesh = std::make_unique<G4DNAMesh>(fpMesh->GetBoundingBox(), fPixel);
228 <<
229 if(!CheckingReactionRadius(fpMesh->GetReso <<
230 { <<
231 G4String WarMessage = "resolution is not <<
232 std::to_string(fpM <<
233 G4Exception("G4DNAEventScheduler::Initia <<
234 JustWarning, WarMessage); <<
235 } <<
236 248
237 // Scavenger(); 249 // Scavenger();
238 250
239 auto pScavengerMaterial = dynamic_cast<G4D 251 auto pScavengerMaterial = dynamic_cast<G4DNAScavengerMaterial*>(
240 G4Scheduler::Instance()->GetScavengerMat 252 G4Scheduler::Instance()->GetScavengerMaterial());
241 if(pScavengerMaterial == nullptr) 253 if(pScavengerMaterial == nullptr)
242 { 254 {
243 G4cout << "There is no scavenger" << G4e << 255 G4cout << "pScavengerMaterial == nullptr" << G4endl;
244 } 256 }
245 else 257 else
246 { 258 {
247 if(fVerbose > 1) 259 if(fVerbose > 1)
248 { 260 {
249 pScavengerMaterial->PrintInfo(); 261 pScavengerMaterial->PrintInfo();
250 } 262 }
251 } 263 }
252 264
253 Voxelizing(); 265 Voxelizing();
254 fpGillespieReaction->SetVoxelMesh(*fpMesh) 266 fpGillespieReaction->SetVoxelMesh(*fpMesh);
255 fpGillespieReaction->SetEventSet(fpEventSe 267 fpGillespieReaction->SetEventSet(fpEventSet.get());
256 fpGillespieReaction->SetTimeStep(0);// res << 268 fpGillespieReaction->SetTimeStep(
>> 269 0); // reset fTimeStep = 0 in fpGillespieReaction
257 fpGillespieReaction->Initialize(); 270 fpGillespieReaction->Initialize();
258 fpGillespieReaction->CreateEvents(); <<
259 fpUpdateSystem->SetMesh(fpMesh.get()); 271 fpUpdateSystem->SetMesh(fpMesh.get());
260 ClearAndReChargeCounter(); 272 ClearAndReChargeCounter();
261 fInitialized = true; 273 fInitialized = true;
262 } 274 }
263 275
264 if(fVerbose > 0) 276 if(fVerbose > 0)
265 { 277 {
266 fpUpdateSystem->SetVerbose(1); 278 fpUpdateSystem->SetVerbose(1);
267 } 279 }
268 280
269 if(fVerbose > 2) 281 if(fVerbose > 2)
270 { 282 {
271 fpMesh->PrintMesh(); 283 fpMesh->PrintMesh();
272 } 284 }
273 } 285 }
274 void G4DNAEventScheduler::InitializeInMesh() 286 void G4DNAEventScheduler::InitializeInMesh()
275 { 287 {
276 if(fPixel <= 1) 288 if(fPixel <= 1)
277 { 289 {
278 fRunning = false; 290 fRunning = false;
279 return; 291 return;
280 } 292 }
281 // TEST /3 293 // TEST /3
282 ReVoxelizing(fPixel / 2); // 294 ReVoxelizing(fPixel / 2); //
283 // ReVoxelizing(fPixel/3);// 295 // ReVoxelizing(fPixel/3);//
284 296
285 fpGillespieReaction->SetVoxelMesh(*fpMesh); 297 fpGillespieReaction->SetVoxelMesh(*fpMesh);
286 fpUpdateSystem->SetMesh(fpMesh.get()); 298 fpUpdateSystem->SetMesh(fpMesh.get());
287 fpGillespieReaction->CreateEvents(); << 299 fpGillespieReaction->Initialize();
288 } 300 }
289 301
290 void G4DNAEventScheduler::ResetInMesh() 302 void G4DNAEventScheduler::ResetInMesh()
291 { 303 {
292 if(fVerbose > 0) 304 if(fVerbose > 0)
293 { 305 {
294 G4cout 306 G4cout
295 << "*** End Processing In Mesh and reset 307 << "*** End Processing In Mesh and reset Mesh, EventSet for new Mesh!!!!"
296 << G4endl; 308 << G4endl;
297 } 309 }
298 fpEventSet->RemoveEventSet(); 310 fpEventSet->RemoveEventSet();
299 fInitialized = false; 311 fInitialized = false;
300 fIsChangeMesh = false; 312 fIsChangeMesh = false;
301 fReactionNumber = 0; 313 fReactionNumber = 0;
302 fJumpingNumber = 0; 314 fJumpingNumber = 0;
303 fStepNumberInMesh = 0; 315 fStepNumberInMesh = 0;
304 } 316 }
305 317
306 G4double G4DNAEventScheduler::GetStartTime() c 318 G4double G4DNAEventScheduler::GetStartTime() const { return fStartTime; }
307 319
308 G4double G4DNAEventScheduler::GetEndTime() con 320 G4double G4DNAEventScheduler::GetEndTime() const { return fEndTime; }
309 321
310 [[maybe_unused]] G4double G4DNAEventScheduler: << 322 [[maybe_unused]] G4double G4DNAEventScheduler::GetTimeStep() const { return fTimeStep; }
311 { <<
312 return fTimeStep; <<
313 } <<
314 323
315 G4int G4DNAEventScheduler::GetVerbose() const 324 G4int G4DNAEventScheduler::GetVerbose() const { return fVerbose; }
316 325
317 [[maybe_unused]] void G4DNAEventScheduler::Set << 326 [[maybe_unused]] void G4DNAEventScheduler::SetMaxNbSteps(G4int max) { fMaxStep = max; }
318 { <<
319 fMaxStep = max; <<
320 } <<
321 327
322 [[maybe_unused]] void G4DNAEventScheduler::Set 328 [[maybe_unused]] void G4DNAEventScheduler::SetStartTime(G4double time)
323 { 329 {
324 fStartTime = time; << 330 fStartTime = time;
325 fGlobalTime = fStartTime; <<
326 } 331 }
327 332
328 void G4DNAEventScheduler::Stop() { fRunning = 333 void G4DNAEventScheduler::Stop() { fRunning = false; }
329 void G4DNAEventScheduler::Run() 334 void G4DNAEventScheduler::Run()
330 { 335 {
331 G4Timer localtimer; 336 G4Timer localtimer;
332 if(fVerbose > 2) << 337 if(fVerbose > 0)
333 { 338 {
334 localtimer.Start(); 339 localtimer.Start();
335 G4cout << "***G4DNAEventScheduler::Run*** 340 G4cout << "***G4DNAEventScheduler::Run*** for Pixel : " << fPixel << G4endl;
336 } 341 }
337 while(fEndTime > fGlobalTime && fRunning) 342 while(fEndTime > fGlobalTime && fRunning)
338 { 343 {
339 RunInMesh(); 344 RunInMesh();
340 } 345 }
341 if(fVerbose > 2) << 346 if(fVerbose > 0)
342 { 347 {
343 if(!fRunning) 348 if(!fRunning)
344 { 349 {
345 G4cout << " StepNumber(" << fStepNumber 350 G4cout << " StepNumber(" << fStepNumber << ") = MaxStep(" << fMaxStep
346 << ")" << G4endl; 351 << ")" << G4endl;
347 } 352 }
348 else if(fEndTime <= fGlobalTime) 353 else if(fEndTime <= fGlobalTime)
349 { 354 {
350 G4cout << " GlobalTime(" << fGlobalTime 355 G4cout << " GlobalTime(" << fGlobalTime << ") > EndTime(" << fEndTime
351 << ")" 356 << ")"
352 << " StepNumber : " << fStepNumbe 357 << " StepNumber : " << fStepNumber << G4endl;
353 } 358 }
354 localtimer.Stop(); 359 localtimer.Stop();
355 G4cout << "***G4DNAEventScheduler::Ending: 360 G4cout << "***G4DNAEventScheduler::Ending::"
356 << G4BestUnit(fGlobalTime, "Time") 361 << G4BestUnit(fGlobalTime, "Time")
357 << " Events left : " << fpEventSet- 362 << " Events left : " << fpEventSet->size() << G4endl;
358 if(fVerbose > 1) << 363 if(fVerbose > 1) {
359 { <<
360 fpMesh->PrintMesh(); 364 fpMesh->PrintMesh();
361 } 365 }
362 G4cout << " Computing Time : " << localtim 366 G4cout << " Computing Time : " << localtimer << G4endl;
363 } 367 }
364 Reset(); 368 Reset();
365 } 369 }
366 370
367 void G4DNAEventScheduler::RunInMesh() 371 void G4DNAEventScheduler::RunInMesh()
368 { 372 {
369 if(!fInitialized) 373 if(!fInitialized)
370 { 374 {
371 InitializeInMesh(); 375 InitializeInMesh();
372 } 376 }
373 if(fVerbose > 0) << 377 G4Timer localtimerInMesh;
>> 378 // if (fVerbose > 0)
374 { 379 {
375 G4double resolution = fpMesh->GetResolutio << 380 localtimerInMesh.Start();
376 G4cout << "At Time : " << std::setw(7) << << 381 G4double C = 20;
377 << " the Mesh has " << fPixel << " << 382 G4double D = G4MoleculeTable::Instance()
378 << " voxels with Resolution " << G4 << 383 ->GetConfiguration("H2O2")
379 << " during next " << 384 ->GetDiffusionCoefficient();
380 << G4BestUnit(resolution * resoluti << 385 G4double transferTime = std::pow(fpMesh->GetResolution(), 2) * C / (6 * D);
>> 386 G4cout << "***G4DNAEventScheduler::RunInMesh*** for Pixel : " << fPixel
>> 387 << " transferTime : " << G4BestUnit(transferTime, "Time") << G4endl;
>> 388 G4cout << " resolution : " << G4BestUnit(fpMesh->GetResolution(), "Length")
381 << G4endl; 389 << G4endl;
382 } 390 }
383 391
384 if(fVerbose > 2) 392 if(fVerbose > 2)
385 { 393 {
386 fpMesh->PrintMesh(); 394 fpMesh->PrintMesh();
387 } 395 }
388 396
389 if(fpUserMeshAction != nullptr) 397 if(fpUserMeshAction != nullptr)
390 { 398 {
391 fpUserMeshAction->BeginOfMesh(fpMesh.get() 399 fpUserMeshAction->BeginOfMesh(fpMesh.get(), fGlobalTime);
392 } 400 }
393 401
394 // if diffusive jumping is avaiable, EventSe 402 // if diffusive jumping is avaiable, EventSet is never empty
395 <<
396 while(!fpEventSet->Empty() && !fIsChangeMesh 403 while(!fpEventSet->Empty() && !fIsChangeMesh && fEndTime > fGlobalTime)
397 { 404 {
398 Stepping(); 405 Stepping();
399 fGlobalTime = fTimeStep + fStartTime; 406 fGlobalTime = fTimeStep + fStartTime;
400 407
401 if(fpUserMeshAction != nullptr) 408 if(fpUserMeshAction != nullptr)
402 { 409 {
403 fpUserMeshAction->InMesh(fpMesh.get(), f 410 fpUserMeshAction->InMesh(fpMesh.get(), fGlobalTime);
404 } 411 }
405 412
406 if(fVerbose > 2) 413 if(fVerbose > 2)
407 { 414 {
408 G4cout << "fGlobalTime : " << G4BestUnit 415 G4cout << "fGlobalTime : " << G4BestUnit(fGlobalTime, "Time")
409 << " fTimeStep : " << G4BestUnit( 416 << " fTimeStep : " << G4BestUnit(fTimeStep, "Time") << G4endl;
410 } 417 }
411 G4double resolution = fpMesh->GetResolutio << 418
412 fTransferTime = resolution * resolut << 419 G4double C = 20;
413 if(fTransferTime == 0) << 420 G4double D = G4MoleculeTable::Instance()
>> 421 ->GetConfiguration("H2O2")
>> 422 ->GetDiffusionCoefficient();
>> 423 if(D == 0)
414 { 424 {
415 G4ExceptionDescription exceptionDescript 425 G4ExceptionDescription exceptionDescription;
416 exceptionDescription << "fTransferTime = << 426 exceptionDescription << "D == 0";
417 G4Exception("G4DNAEventScheduler::RunInM 427 G4Exception("G4DNAEventScheduler::RunInMesh", "G4DNAEventScheduler001",
418 FatalErrorInArgument, except 428 FatalErrorInArgument, exceptionDescription);
419 } 429 }
420 if(fTransferTime < fTimeStep && << 430 G4double transferTime = std::pow(fpMesh->GetResolution(), 2) * C / (6 * D);
421 fPixel != 1) // dont change Mesh if fP << 431
>> 432 // if(fStepNumberInMesh > 40000 && fPixel != 1)
>> 433 if(transferTime < fTimeStep &&
>> 434 fPixel != 1) // dont change Mesj if fPixel == 1
422 { 435 {
>> 436 if(fVerbose > 1)
>> 437 {
>> 438 G4cout << " Pixels : " << fPixel << " resolution : "
>> 439 << G4BestUnit(fpMesh->GetResolution(), "Length")
>> 440 << " fStepNumberInMesh : " << fStepNumberInMesh
>> 441 << " at fGlobalTime : " << G4BestUnit(fGlobalTime, "Time")
>> 442 << " at fTimeStep : " << G4BestUnit(fTimeStep, "Time")
>> 443 << " fReactionNumber : " << fReactionNumber
>> 444 << " transferTime : " << G4BestUnit(transferTime, "Time")
>> 445 << G4endl;
>> 446 }
423 if(fSetChangeMesh) 447 if(fSetChangeMesh)
424 { 448 {
425 if(fVerbose > 1) <<
426 { <<
427 G4cout << " Pixels : " << fPixel << <<
428 << G4BestUnit(fpMesh->GetReso <<
429 << " fStepNumberInMesh : " < <<
430 << " at fGlobalTime : " << G4 <<
431 << " at fTimeStep : " << G4Be <<
432 << " fReactionNumber : " << <<
433 << " transferTime : " << G4Be <<
434 << G4endl; <<
435 } <<
436 fIsChangeMesh = true; 449 fIsChangeMesh = true;
437 } 450 }
438 } 451 }
439 } 452 }
440 453
441 if(fVerbose > 1) 454 if(fVerbose > 1)
442 { 455 {
>> 456 localtimerInMesh.Stop();
443 G4cout << "***G4DNAEventScheduler::Ending: 457 G4cout << "***G4DNAEventScheduler::Ending::"
444 << G4BestUnit(fGlobalTime, "Time") 458 << G4BestUnit(fGlobalTime, "Time")
445 << " Event left : " << fpEventSet-> 459 << " Event left : " << fpEventSet->size() << G4endl;
446 G4cout << " Due to : "; << 460 G4cout << " Computing Time : " << localtimerInMesh << " Due to : ";
447 if(fpEventSet->Empty()) 461 if(fpEventSet->Empty())
448 { 462 {
449 G4cout << "EventSet is Empty" << G4endl; 463 G4cout << "EventSet is Empty" << G4endl;
450 } 464 }
451 else if(fIsChangeMesh) 465 else if(fIsChangeMesh)
452 { 466 {
453 G4cout << "Changing Mesh from : " << fPi 467 G4cout << "Changing Mesh from : " << fPixel
454 << " pixels to : " << fPixel / 2 468 << " pixels to : " << fPixel / 2 << " pixels" << G4endl;
455 G4cout << "Info : ReactionNumber : " << 469 G4cout << "Info : ReactionNumber : " << fReactionNumber
456 << " JumpingNumber : " << fJump 470 << " JumpingNumber : " << fJumpingNumber << G4endl;
457 } 471 }
458 else if(fEndTime > fGlobalTime) 472 else if(fEndTime > fGlobalTime)
459 { 473 {
460 G4cout << " GlobalTime(" << fGlobalTime 474 G4cout << " GlobalTime(" << fGlobalTime << ") > EndTime(" << fEndTime
461 << ")" 475 << ")"
462 << " StepNumber : " << fStepNumbe 476 << " StepNumber : " << fStepNumber << G4endl;
463 } 477 }
464 if(fVerbose > 2) 478 if(fVerbose > 2)
465 { 479 {
466 fpMesh->PrintMesh(); 480 fpMesh->PrintMesh();
467 } 481 }
468 G4cout << G4endl; 482 G4cout << G4endl;
469 } 483 }
470 484
471 if(fpUserMeshAction != nullptr) 485 if(fpUserMeshAction != nullptr)
472 { 486 {
473 fpUserMeshAction->EndOfMesh(fpMesh.get(), 487 fpUserMeshAction->EndOfMesh(fpMesh.get(), fGlobalTime);
474 } 488 }
475 ResetInMesh(); 489 ResetInMesh();
476 } 490 }
477 491
478 void G4DNAEventScheduler::Stepping() // this 492 void G4DNAEventScheduler::Stepping() // this event loop
479 { 493 {
480 fStepNumber < fMaxStep ? fStepNumber++ : sta << 494 fStepNumber < fMaxStep ? fStepNumber++ : fRunning = false;
>> 495
481 if(fpEventSet->size() > fpMesh->size()) 496 if(fpEventSet->size() > fpMesh->size())
482 { 497 {
483 G4ExceptionDescription exceptionDescriptio 498 G4ExceptionDescription exceptionDescription;
484 exceptionDescription << 499 exceptionDescription << "fpEventSet->size() > fpMesh->size()";
485 << "impossible that fpEventSet->size() > <<
486 G4Exception("G4DNAEventScheduler::Stepping 500 G4Exception("G4DNAEventScheduler::Stepping", "G4DNAEventScheduler002",
487 FatalErrorInArgument, exceptio 501 FatalErrorInArgument, exceptionDescription);
488 } << 502 };
489 503
490 auto selected = fpEventSet->begin(); << 504 auto selected = fpEventSet->begin();
491 auto index = (*selected)->GetIndex(); << 505 const auto& key = (*selected)->GetKey();
>> 506 auto index = fpMesh->GetIndex(key);
492 507
493 if(fVerbose > 1) 508 if(fVerbose > 1)
494 { 509 {
495 G4cout << "G4DNAEventScheduler::Stepping() 510 G4cout << "G4DNAEventScheduler::Stepping()*********************************"
496 "*******" 511 "*******"
497 << G4endl; 512 << G4endl;
498 (*selected)->PrintEvent(); 513 (*selected)->PrintEvent();
499 } 514 }
500 515
501 // get selected time step 516 // get selected time step
502 fTimeStep = (*selected)->GetTime(); 517 fTimeStep = (*selected)->GetTime();
503 518
504 // selected data 519 // selected data
505 auto pJumping = (*selected)->GetJumpingData 520 auto pJumping = (*selected)->GetJumpingData();
506 auto pReaction = (*selected)->GetReactionDat 521 auto pReaction = (*selected)->GetReactionData();
507 522
508 fpUpdateSystem->SetGlobalTime(fTimeStep + 523 fpUpdateSystem->SetGlobalTime(fTimeStep +
509 fStartTime); 524 fStartTime); // this is just for printing
>> 525
>> 526 if(pJumping == nullptr && pReaction == nullptr)
>> 527 {
>> 528 G4ExceptionDescription exceptionDescription;
>> 529 exceptionDescription << "pJumping == nullptr && pReaction == nullptr";
>> 530 G4Exception("G4DNAEventScheduler::Stepping", "G4DNAEventScheduler003",
>> 531 FatalErrorInArgument, exceptionDescription);
>> 532 }
>> 533
510 fpGillespieReaction->SetTimeStep(fTimeStep); 534 fpGillespieReaction->SetTimeStep(fTimeStep);
511 if(pJumping == nullptr && pReaction != nullp << 535
>> 536 if(pJumping == nullptr)
512 { 537 {
513 fpUpdateSystem->UpdateSystem(index, *pReac 538 fpUpdateSystem->UpdateSystem(index, *pReaction);
514 fpEventSet->RemoveEvent(selected); <<
515 <<
516 //hoang's exp <<
517 if(fpGillespieReaction->SetEquilibrium(pRe <<
518 { <<
519 ResetEventSet(); <<
520 } <<
521 //end Hoang's exp <<
522 539
>> 540 fpEventSet->RemoveEvent(selected);
523 // create new event 541 // create new event
524 fpGillespieReaction->CreateEvent(index); << 542 fpGillespieReaction->CreateEvent(key);
525 fReactionNumber++; 543 fReactionNumber++;
>> 544
526 // recordTime in reaction 545 // recordTime in reaction
527 RecordTime(); 546 RecordTime();
528 } 547 }
529 else if(pJumping != nullptr && pReaction == << 548 else if(pReaction == nullptr)
530 { 549 {
531 // dont change this 550 // dont change this
532 fpUpdateSystem->UpdateSystem(index, *pJump 551 fpUpdateSystem->UpdateSystem(index, *pJumping);
533 // save jumping Index before delete select << 552 auto jumpingKey = fpMesh->GetKey(pJumping->second);
534 auto jumpingIndex = pJumping->second; <<
535 fpEventSet->RemoveEvent(selected); 553 fpEventSet->RemoveEvent(selected);
>> 554
536 // create new event 555 // create new event
537 // should create Jumping before key 556 // should create Jumping before key
538 fpGillespieReaction->CreateEvent(jumpingIn << 557 fpGillespieReaction->CreateEvent(jumpingKey);
539 fpGillespieReaction->CreateEvent(index); << 558 fpGillespieReaction->CreateEvent(key);
>> 559
540 fJumpingNumber++; 560 fJumpingNumber++;
541 } 561 }
542 else <<
543 { <<
544 G4ExceptionDescription exceptionDescriptio <<
545 exceptionDescription << "pJumping == nullp <<
546 G4Exception("G4DNAEventScheduler::Stepping <<
547 FatalErrorInArgument, exceptio <<
548 } <<
549 <<
550 if(fVerbose > 1) 562 if(fVerbose > 1)
551 { 563 {
552 G4cout << "G4DNAEventScheduler::Stepping:: 564 G4cout << "G4DNAEventScheduler::Stepping::end "
553 "Print************************** 565 "Print***********************************"
554 << G4endl; 566 << G4endl;
555 G4cout << G4endl; 567 G4cout << G4endl;
556 } 568 }
557 fStepNumberInMesh++; 569 fStepNumberInMesh++;
558 } 570 }
559 571
560 void G4DNAEventScheduler::SetEndTime(const G4d 572 void G4DNAEventScheduler::SetEndTime(const G4double& endTime)
561 { 573 {
562 fEndTime = endTime; 574 fEndTime = endTime;
563 } 575 }
564 576
565 void G4DNAEventScheduler::RecordTime() 577 void G4DNAEventScheduler::RecordTime()
566 { 578 {
567 auto recordTime = *fLastRecoredTime; 579 auto recordTime = *fLastRecoredTime;
568 if(fGlobalTime >= recordTime && fCounterMap[ 580 if(fGlobalTime >= recordTime && fCounterMap[recordTime].empty())
569 { 581 {
570 auto begin = fpMesh->begin(); 582 auto begin = fpMesh->begin();
571 auto end = fpMesh->end(); 583 auto end = fpMesh->end();
572 for(; begin != end; begin++) 584 for(; begin != end; begin++)
573 { 585 {
574 const auto& mapData = std::get<2>(*begin << 586 auto node = begin->second;
575 if(mapData.empty()) << 587 if(node == nullptr) {
576 { <<
577 continue; 588 continue;
578 } 589 }
579 for(const auto& it : mapData) << 590 for(const auto& it : node->GetMapList())
580 { 591 {
581 fCounterMap[recordTime][it.first] += i 592 fCounterMap[recordTime][it.first] += it.second;
582 } 593 }
583 } 594 }
584 fLastRecoredTime++; 595 fLastRecoredTime++;
>> 596
>> 597 #ifdef DEBUG
>> 598 PrintRecordTime();
>> 599 G4MoleculeTable* pMoleculeTable = G4MoleculeTable::Instance();
>> 600 auto iter = pMoleculeTable->GetConfigurationIterator();
>> 601 iter.reset();
>> 602 while(iter())
>> 603 {
>> 604 auto conf = iter.value();
>> 605
>> 606 G4cout << "GlobalTime : " << G4BestUnit(fGlobalTime, "Time")
>> 607 << " recordTime : " << G4BestUnit(recordTime, "Time") << " "
>> 608 << conf->GetName()
>> 609 << " number : " << fCounterMap[recordTime][conf]
>> 610 << " MoleculeCounter : "
>> 611 << G4MoleculeCounter::Instance()->GetCurrentNumberOf(conf)
>> 612 << G4endl;
>> 613
>> 614 assert(G4MoleculeCounter::Instance()->GetCurrentNumberOf(conf) ==
>> 615 fCounterMap[recordTime][conf]);
>> 616 }
>> 617 #endif
585 } 618 }
586 } 619 }
587 620
588 void G4DNAEventScheduler::PrintRecordTime() 621 void G4DNAEventScheduler::PrintRecordTime()
589 { 622 {
590 G4cout << "CounterMap.size : " << fCounterMa << 623 G4cout << "fCounterMap.size : " << fCounterMap.size() << G4endl;
>> 624
591 for(const auto& i : fCounterMap) 625 for(const auto& i : fCounterMap)
592 { 626 {
593 auto map = i.second; 627 auto map = i.second;
594 auto begin = map.begin(); // 628 auto begin = map.begin(); //
595 auto end = map.end(); // 629 auto end = map.end(); //
596 for(; begin != end; begin++) 630 for(; begin != end; begin++)
597 { 631 {
598 auto molecule = begin->first; 632 auto molecule = begin->first;
599 auto number = begin->second; 633 auto number = begin->second;
600 if(number == 0) 634 if(number == 0)
601 { 635 {
602 continue; 636 continue;
603 } 637 }
604 G4cout << "molecule : " << molecule->Get 638 G4cout << "molecule : " << molecule->GetName() << " number : " << number
605 << G4endl; 639 << G4endl;
606 } 640 }
607 } 641 }
608 } 642 }
609 643
610 std::map<G4double /*time*/, G4DNAEventSchedule 644 std::map<G4double /*time*/, G4DNAEventScheduler::MapCounter>
611 G4DNAEventScheduler::GetCounterMap() const 645 G4DNAEventScheduler::GetCounterMap() const
612 { 646 {
613 return fCounterMap; 647 return fCounterMap;
614 } 648 }
615 649
616 void G4DNAEventScheduler::SetUserMeshAction( 650 void G4DNAEventScheduler::SetUserMeshAction(
617 std::unique_ptr<G4UserMeshAction> pUserMeshA 651 std::unique_ptr<G4UserMeshAction> pUserMeshAction)
618 { 652 {
619 fpUserMeshAction = std::move(pUserMeshAction 653 fpUserMeshAction = std::move(pUserMeshAction);
620 } 654 }
621 655
622 G4DNAMesh* G4DNAEventScheduler::GetMesh() cons 656 G4DNAMesh* G4DNAEventScheduler::GetMesh() const { return fpMesh.get(); }
623 657
624 G4int G4DNAEventScheduler::GetPixels() const { 658 G4int G4DNAEventScheduler::GetPixels() const { return fPixel; }
625 659
626 G4bool G4DNAEventScheduler::CheckingReactionRa 660 G4bool G4DNAEventScheduler::CheckingReactionRadius(G4double resolution)
627 { 661 {
628 auto pMolecularReactionTable = G4DNAMolecula 662 auto pMolecularReactionTable = G4DNAMolecularReactionTable::Instance();
629 auto reactionDataList = pMolecularReactionTa 663 auto reactionDataList = pMolecularReactionTable->GetVectorOfReactionData();
630 if(reactionDataList.empty()) 664 if(reactionDataList.empty())
631 { 665 {
632 G4cout << "reactionDataList.empty()" << G4 666 G4cout << "reactionDataList.empty()" << G4endl;
633 return true; 667 return true;
634 } 668 }
635 << 669 else
636 for(auto it : reactionDataList) <<
637 { <<
638 if(it->GetEffectiveReactionRadius() >= res <<
639 { <<
640 G4cout << it->GetReactant1()->GetName() <<
641 << it->GetReactant2()->GetName() <<
642 G4cout << "G4DNAEventScheduler::Reaction <<
643 << G4BestUnit(it->GetEffectiveRea <<
644 << G4endl; <<
645 G4cout << "resolution : " << G4BestUnit( <<
646 return false; <<
647 } <<
648 } <<
649 return true; <<
650 } <<
651 <<
652 void G4DNAEventScheduler::ResetEventSet() <<
653 { <<
654 fpEventSet->RemoveEventSet(); <<
655 fpGillespieReaction->CreateEvents(); <<
656 } <<
657 <<
658 void G4DNAEventScheduler::LastRegisterForCount <<
659 { <<
660 if(fLastRecoredTime == fTimeToRecord.end()) <<
661 { <<
662 //fully recorded -> happy ending <<
663 return; <<
664 }else <<
665 { 670 {
666 auto lastRecorded = --fLastRecoredTime;//f << 671 for(auto it : reactionDataList)
667 while (fLastRecoredTime != fTimeToRecord.e <<
668 { 672 {
669 fCounterMap[*fLastRecoredTime] = fCounte << 673 if(it->GetEffectiveReactionRadius() >= resolution / CLHEP::pi)
670 fLastRecoredTime++; << 674 {
>> 675 G4cout << it->GetReactant1()->GetName() << " + "
>> 676 << it->GetReactant2()->GetName() << G4endl;
>> 677 G4cout << "G4DNAEventScheduler::ReactionRadius : "
>> 678 << G4BestUnit(it->GetEffectiveReactionRadius(), "Length")
>> 679 << G4endl;
>> 680 G4cout << "resolution : " << G4BestUnit(resolution, "Length") << G4endl;
>> 681 return false;
>> 682 }
671 } 683 }
>> 684 return true;
672 } 685 }
673 <<
674 } 686 }