| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25
26 #include "Scorer.hh"
27
28 #include "PrimaryGeneratorAction.hh"
29 #include "TimeStepAction.hh"
30
31 #include "G4AnalysisManager.hh"
32 #include "G4DNAEventScheduler.hh"
33 #include "G4DNAScavengerMaterial.hh"
34 #include "G4Event.hh"
35 #include "G4MoleculeTable.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4RunManager.hh"
38 #include "G4Scheduler.hh"
39 #include "G4UIcmdWithADoubleAndUnit.hh"
40 #include "G4UIcmdWithAnInteger.hh"
41 #include "G4VChemistryWorld.hh"
42
43 #include <G4EventManager.hh>
44 #include <G4MolecularConfiguration.hh>
45 #include <G4MoleculeCounter.hh>
46 #include <G4SystemOfUnits.hh>
47 #include <globals.hh>
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49
50 // Dose
51 Dose::Dose()
52 : G4UImessenger(),
53 fpDoseDir(new G4UIdirectory("/scorer/Dose/")),
54 fpAddDoseCutOff(new G4UIcmdWithADoubleAndUnit("/scorer/Dose/cutoff", this)),
55 fpAddDoseToAbort(new G4UIcmdWithADoubleAndUnit("/scorer/Dose/abortedDose", this))
56 {
57 fpDoseDir->SetGuidance("Dose scorer commands");
58 }
59
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61
62 void Dose::SetNewValue(G4UIcommand* command, G4String newValue)
63 {
64 if (command == fpAddDoseCutOff.get()) {
65 fDosesCutOff = G4UIcmdWithADoubleAndUnit::GetNewDoubleValue(newValue);
66 }
67 if (command == fpAddDoseToAbort.get()) {
68 fDosesToAbort = G4UIcmdWithADoubleAndUnit::GetNewDoubleValue(newValue);
69 }
70 }
71
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73
74 template<>
75 void Scorer<Dose>::SetChemistryWorld(G4VChemistryWorld* pChemistryWorld)
76 {
77 fpChemistryWorld = pChemistryWorld;
78 }
79
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81
82 template<>
83 G4VChemistryWorld* Scorer<Dose>::GetChemistryWorld() const
84 {
85 return fpChemistryWorld;
86 }
87
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
89
90 template<>
91 void Scorer<Dose>::clear()
92 {
93 fpScorer->fCumulatedDose = 0.;
94 }
95
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
97
98 template<>
99 void Scorer<Dose>::Initialize(G4HCofThisEvent* HCE)
100 {
101 clear();
102 fpEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
103 if (fHCID < 0) {
104 fHCID = GetCollectionID(0);
105 }
106 HCE->AddHitsCollection(fHCID, (G4VHitsCollection*)fpEvtMap);
107 }
108
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
110
111 template<>
112 void Scorer<Dose>::EndOfEvent(G4HCofThisEvent*)
113 {
114 if (!G4RunManager::GetRunManager()->GetCurrentEvent()->IsAborted()) {
115 fpEvtMap->add(0, fpScorer->fDosesCutOff);
116 }
117 fpScorer->fCumulatedDose = 0.;
118 }
119
120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
121
122 template<>
123 G4bool Scorer<Dose>::ProcessHits(G4Step* aStep, G4TouchableHistory*)
124 {
125 auto currentEvent = G4EventManager::GetEventManager();
126 const G4Track* track = aStep->GetTrack();
127 auto boundingBox = fpChemistryWorld->GetChemistryBoundary();
128 G4double V = boundingBox->Volume() / cm3;
129 G4double edep = aStep->GetTotalEnergyDeposit();
130 if (edep == 0.) {
131 return false;
132 }
133 (fpScorer->fCumulatedDose) += edep;
134 if (track->GetParentID() == 0 && aStep->IsFirstStepInVolume()) {
135 G4double DoseInGray = ((fpScorer->fCumulatedDose) / eV) / (0.001 * V * 6.242e+18);
136 if (DoseInGray > fpScorer->fDosesCutOff / gray) {
137 G4cout << "_____________________________________________________________________________"
138 << G4endl;
139 auto name = currentEvent->GetConstCurrentEvent()
140 ->GetPrimaryVertex()
141 ->GetPrimary()
142 ->GetParticleDefinition()
143 ->GetParticleName();
144 auto energy =
145 currentEvent->GetConstCurrentEvent()->GetPrimaryVertex()->GetPrimary()->GetKineticEnergy();
146 G4cout << "Stop this beam line (" << name << ", " << energy
147 << " MeV) at actual dose: " << DoseInGray
148 << " Gy. Cut-off dose: " << fpScorer->fDosesCutOff / gray << " Gy" << G4endl;
149 G4cout << "The beam of " << 1000000 - currentEvent->GetStackManager()->GetNUrgentTrack() - 1
150 << " tracks"
151 << " in a volume of " << V * 1e+12 // convert cm3 to um3
152 << " um3. Total deposit energy: " << fpScorer->fCumulatedDose / eV << " eV. "
153 << G4endl;
154 if (DoseInGray > fpScorer->fDosesToAbort / gray) {
155 G4cout << "Abort this beam line (" << name << ", " << energy
156 << " MeV) at actual dose: " << DoseInGray << " Gy." << G4endl;
157 G4RunManager::GetRunManager()->AbortEvent();
158 }
159 G4cout << "_____________________________________________________________________________"
160 << G4endl;
161 auto myTrack = ((G4Track*)track);
162 myTrack->SetTrackStatus(fStopAndKill);
163 auto secondaries = track->GetStep()->GetSecondaryInCurrentStep();
164 if (!secondaries->empty()) {
165 for (auto it : *(secondaries)) {
166 if (it != nullptr) {
167 ((G4Track*)it)->SetTrackStatus(fStopAndKill);
168 }
169 }
170 }
171 currentEvent->GetStackManager()->ClearUrgentStack();
172 }
173 }
174 return true;
175 }
176
177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
178
179 // Gvalues
180 Gvalues::Gvalues()
181 : G4UImessenger(),
182 fTimeLimit(G4Scheduler::Instance()->GetEndTime()),
183 fSpeciesdir(new G4UIdirectory("/scorer/Gvalues/")),
184 fTimeBincmd(new G4UIcmdWithAnInteger("/scorer/Gvalues/nOfTimeBins", this)),
185 fAddTimeToRecordcmd(new G4UIcmdWithADoubleAndUnit("/scorer/Gvalues/addTimeToRecord", this))
186 {
187 fSpeciesdir->SetGuidance("ScoreSpecies commands");
188 G4MoleculeCounter::Instance()->ResetCounter();
189 }
190
191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
192
193 void Gvalues::SetNewValue(G4UIcommand* command, G4String newValue)
194 {
195 if (command == fAddTimeToRecordcmd.get()) {
196 G4double cmdTime = fAddTimeToRecordcmd->GetNewDoubleValue(newValue);
197 if (fTimeLimit >= cmdTime) {
198 AddTimeToRecord(cmdTime);
199 }
200 else {
201 AddTimeToRecord(fTimeLimit);
202 }
203 }
204 if (command == fTimeBincmd.get()) {
205 G4int cmdBins = fTimeBincmd->GetNewIntValue(newValue);
206 G4double timeMin = 1 * ps;
207 G4double timeMax = G4Scheduler::Instance()->GetEndTime() - timeMin;
208 G4double timeLogMin = std::log10(timeMin);
209 G4double timeLogMax = std::log10(timeMax);
210 for (int i = 0; i <= cmdBins; i++) {
211 AddTimeToRecord(std::pow(10, timeLogMin + i * (timeLogMax - timeLogMin) / cmdBins));
212 }
213 }
214 }
215
216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
217
218 void Gvalues::WriteWithAnalysisManager(G4VAnalysisManager* analysisManager, const std::string& out)
219 {
220 analysisManager->CreateNtuple(out, out);
221 G4cout << "NtupleID : " << fRunID << " name : " << out << G4endl;
222 analysisManager->CreateNtupleIColumn(fRunID, "speciesID");
223 analysisManager->CreateNtupleIColumn(fRunID, "number");
224 analysisManager->CreateNtupleIColumn(fRunID, "nEvent");
225 analysisManager->CreateNtupleSColumn(fRunID, "speciesName");
226 analysisManager->CreateNtupleDColumn(fRunID, "time");
227 analysisManager->CreateNtupleDColumn(fRunID, "sumG");
228 analysisManager->CreateNtupleDColumn(fRunID, "sumG2");
229 analysisManager->FinishNtuple(fRunID);
230
231 for (const auto& it_map1 : fSpeciesInfoPerTime) {
232 const InnerSpeciesMap& map2 = it_map1.second;
233 for (auto it_map2 : map2) {
234 double time = it_map1.first;
235 auto species = it_map2.first;
236 const G4String& name = species->GetName();
237 int molID = it_map2.first->GetMoleculeID();
238 auto number = it_map2.second.fNumber;
239 double G = it_map2.second.fG;
240 double G2 = it_map2.second.fG2;
241
242 analysisManager->FillNtupleIColumn(fRunID, 0, molID); // MolID
243 analysisManager->FillNtupleIColumn(fRunID, 1, number); // Number
244 analysisManager->FillNtupleIColumn(fRunID, 2, fNEvent); // Total nb events
245 analysisManager->FillNtupleSColumn(fRunID, 3, name); // molName
246 analysisManager->FillNtupleDColumn(fRunID, 4, time); // time
247 analysisManager->FillNtupleDColumn(fRunID, 5, G); // G
248 analysisManager->FillNtupleDColumn(fRunID, 6, G2); // G2
249 analysisManager->AddNtupleRow(fRunID);
250 }
251 }
252 fRunID++;
253 }
254
255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
256
257 template<>
258 void Scorer<Gvalues>::clear()
259 {
260 fpEvtMap->clear();
261 fpScorer->fNEvent = 0;
262 fpScorer->fEdep = 0;
263 fpScorer->fSpeciesInfoPerTime.clear();
264 }
265
266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
267
268 template<>
269 void Scorer<Gvalues>::Initialize(G4HCofThisEvent* HCE)
270 {
271 fpEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
272 if (fHCID < 0) {
273 fHCID = GetCollectionID(0);
274 }
275 HCE->AddHitsCollection(fHCID, (G4VHitsCollection*)fpEvtMap);
276 }
277
278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
279
280 template<>
281 void Scorer<Gvalues>::SetChemistryWorld(G4VChemistryWorld* pChemistryWorld)
282 {
283 fpChemistryWorld = pChemistryWorld;
284 }
285
286 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
287
288 template<>
289 G4VChemistryWorld* Scorer<Gvalues>::GetChemistryWorld() const
290 {
291 return fpChemistryWorld;
292 }
293
294 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
295
296 template<>
297 G4bool Scorer<Gvalues>::ProcessHits(G4Step* aStep, G4TouchableHistory*)
298 {
299 G4double edep = aStep->GetTotalEnergyDeposit();
300 if (edep == 0.) {
301 return FALSE;
302 }
303 edep *= aStep->GetPreStepPoint()->GetWeight(); // (Particle Weight)
304 G4int index = GetIndex(aStep);
305 fpEvtMap->add(index, edep);
306 (fpScorer->fEdep) += edep;
307 return TRUE;
308 }
309
310 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
311
312 template<>
313 void Scorer<Gvalues>::SaveScavengerChange()
314 {
315 auto pScavengerMaterial =
316 dynamic_cast<G4DNAScavengerMaterial*>(G4Scheduler::Instance()->GetScavengerMaterial());
317 if (pScavengerMaterial == nullptr) {
318 G4ExceptionDescription errMsg;
319 errMsg << "pScavengerMaterial == nullptr";
320 G4Exception("Scorer<Gvalues>::SaveScavengerChange()", "SaveScavengerChange",
321 FatalErrorInArgument, errMsg);
322 }
323 auto scavengerList = pScavengerMaterial->GetScavengerList();
324 auto V = fpChemistryWorld->GetChemistryBoundary()->Volume();
325
326 for (const auto& it : scavengerList) {
327 if (it == G4MoleculeTable::Instance()->GetConfiguration("H2O")
328 || G4MoleculeTable::Instance()->GetConfiguration("H3Op(B)") == it
329 || G4MoleculeTable::Instance()->GetConfiguration("OHm(B)") == it)
330 {
331 continue;
332 }
333 for (auto time_mol : fpScorer->fTimeToRecord) {
334 int64_t n_mol = pScavengerMaterial->GetNMoleculesAtTime(it, time_mol);
335 if (n_mol < 0) {
336 G4ExceptionDescription errMsg;
337 errMsg << "SaveScavengerChange()::N molecules not valid < 0 : " << it->GetName()
338 << " N : " << n_mol << G4endl;
339 G4Exception("", "N<0", FatalException, errMsg);
340 }
341
342 Gvalues::SpeciesInfo& molInfo = fpScorer->fSpeciesInfoPerTime[time_mol][it];
343 molInfo.fNumber += n_mol;
344 G4double gValue = n_mol / (Avogadro * V * 1.0e-6 /*mm3 to L*/);
345 molInfo.fG += gValue;
346 molInfo.fG2 += gValue * gValue;
347 }
348 }
349 }
350
351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
352
353 template<>
354 void Scorer<Gvalues>::SaveMoleculeCounter()
355 {
356 if (fpEventScheduler == nullptr) {
357 G4Exception("fpEventScheduler == nullptr", "Scorer<Gvalues>::SaveMoleculeCounter()",
358 FatalException, "fpEventScheduler == nullptr");
359 }
360 else {
361 auto counterMap = fpEventScheduler->GetCounterMap();
362 if (counterMap.empty()) {
363 if (!G4MoleculeCounter::Instance()->InUse()) {
364 G4Exception("No counter", "Scorer<Gvalues>::SaveMoleculeCounter()", JustWarning,
365 "G4MoleculeCounter::Instance() is not used");
366 return;
367 }
368
369 G4MoleculeCounter::RecordedMolecules species;
370 species = G4MoleculeCounter::Instance()->GetRecordedMolecules();
371 if (species.get() == nullptr) {
372 return;
373 }
374 else if (species->empty()) {
375 G4cout << "No molecule recorded, energy deposited" << G4endl;
376 ++(fpScorer->fNEvent);
377 fpScorer->fEdep = 0.;
378 G4MoleculeCounter::Instance()->ResetCounter();
379 return;
380 }
381 for (auto molecule : *species) {
382 if (molecule == G4MoleculeTable::Instance()->GetConfiguration("O2")) {
383 continue;
384 }
385 for (auto time_mol : fpScorer->fTimeToRecord) {
386 int n_mol = G4MoleculeCounter::Instance()->GetNMoleculesAtTime(molecule, time_mol);
387
388 if (n_mol < 0) {
389 G4ExceptionDescription errMsg;
390 errMsg << "N molecules not valid < 0 " << G4endl;
391 G4Exception("", "N<0", FatalException, errMsg);
392 }
393
394 Gvalues::SpeciesInfo& molInfo = fpScorer->fSpeciesInfoPerTime[time_mol][molecule];
395 molInfo.fNumber += n_mol;
396 G4double gValue = (n_mol / (fpScorer->fEdep / eV)) * 100.;
397 molInfo.fG += gValue;
398 molInfo.fG2 += gValue * gValue;
399 }
400 }
401 }
402 else {
403 for (const auto& map_mol : counterMap) {
404 auto time_mol = map_mol.first;
405 for (auto it_mol : map_mol.second) {
406 auto molecule = it_mol.first;
407 if (molecule == G4MoleculeTable::Instance()->GetConfiguration("O2")) {
408 continue;
409 }
410 int n_mol = it_mol.second;
411
412 if (n_mol < 0) {
413 G4ExceptionDescription errMsg;
414 errMsg << "N molecules not valid < 0 "
415 << " molecule : " << it_mol.first->GetName() << " N : " << n_mol << G4endl;
416 G4Exception("", "N<0", FatalException, errMsg);
417 }
418
419 Gvalues::SpeciesInfo& molInfo = fpScorer->fSpeciesInfoPerTime[time_mol][molecule];
420 molInfo.fNumber += n_mol;
421 G4double gValue = (n_mol / (fpScorer->fEdep / eV)) * 100.;
422 molInfo.fG += gValue;
423 molInfo.fG2 += gValue * gValue;
424 }
425 }
426 }
427 }
428 }
429
430 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
431
432 template<>
433 void Scorer<Gvalues>::EndOfEvent(G4HCofThisEvent*)
434 {
435 if (G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted()) {
436 fpScorer->fEdep = 0.;
437 G4MoleculeCounter::Instance()->ResetCounter();
438 return;
439 }
440
441 SaveScavengerChange();
442 SaveMoleculeCounter();
443
444 ++(fpScorer->fNEvent);
445 fpScorer->fEdep = 0.;
446
447 G4MoleculeCounter::Instance()->ResetCounter();
448 G4MoleculeCounter::Instance()->Use(true);
449 }
450
451 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
452
453 template<>
454 void Scorer<Gvalues>::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
455 {
456 auto right = dynamic_cast<Scorer<Gvalues>*>(dynamic_cast<G4VPrimitiveScorer*>(workerScorer));
457
458 if (right == nullptr) {
459 return;
460 }
461 if (right == this) {
462 return;
463 }
464 auto it_map1 = right->fpScorer->fSpeciesInfoPerTime.begin();
465 auto end_map1 = right->fpScorer->fSpeciesInfoPerTime.end();
466
467 for (; it_map1 != end_map1; ++it_map1) {
468 Gvalues::InnerSpeciesMap& map2 = it_map1->second;
469 auto it_map2 = map2.begin();
470 auto end_map2 = map2.end();
471
472 for (; it_map2 != end_map2; ++it_map2) {
473 Gvalues::SpeciesInfo& molInfo = fpScorer->fSpeciesInfoPerTime[it_map1->first][it_map2->first];
474 molInfo.fNumber += it_map2->second.fNumber;
475 molInfo.fG += it_map2->second.fG;
476 molInfo.fG2 += it_map2->second.fG2;
477 }
478 }
479 right->fpScorer->fSpeciesInfoPerTime.clear();
480 fpScorer->fNEvent += right->fpScorer->fNEvent;
481 right->fpScorer->fNEvent = 0;
482 right->fpScorer->fEdep = 0.;
483 }
484
485 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
486
487 template<>
488 void Scorer<Gvalues>::OutputAndClear(const std::string& dose)
489 {
490 if (G4Threading::IsWorkerThread()) {
491 return;
492 }
493 G4VAnalysisManager* analysisManager = G4AnalysisManager::Instance();
494 if (analysisManager != nullptr) {
495 this->fpScorer->WriteWithAnalysisManager(analysisManager, dose);
496 }
497 fpScorer->fNEvent = 0;
498 fpScorer->fSpeciesInfoPerTime.clear();
499 }
500
501 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
502