File indexing completed on 2025-02-23 09:22:12
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026 #include "Scorer.hh"
0027
0028 #include "PrimaryGeneratorAction.hh"
0029 #include "TimeStepAction.hh"
0030
0031 #include "G4AnalysisManager.hh"
0032 #include "G4DNAEventScheduler.hh"
0033 #include "G4DNAScavengerMaterial.hh"
0034 #include "G4Event.hh"
0035 #include "G4MoleculeTable.hh"
0036 #include "G4PhysicalConstants.hh"
0037 #include "G4RunManager.hh"
0038 #include "G4Scheduler.hh"
0039 #include "G4UIcmdWithADoubleAndUnit.hh"
0040 #include "G4UIcmdWithAnInteger.hh"
0041 #include "G4VChemistryWorld.hh"
0042
0043 #include <G4EventManager.hh>
0044 #include <G4MolecularConfiguration.hh>
0045 #include <G4MoleculeCounter.hh>
0046 #include <G4SystemOfUnits.hh>
0047 #include <globals.hh>
0048
0049
0050
0051 Dose::Dose()
0052 : G4UImessenger(),
0053 fpDoseDir(new G4UIdirectory("/scorer/Dose/")),
0054 fpAddDoseCutOff(new G4UIcmdWithADoubleAndUnit("/scorer/Dose/cutoff", this)),
0055 fpAddDoseToAbort(new G4UIcmdWithADoubleAndUnit("/scorer/Dose/abortedDose", this))
0056 {
0057 fpDoseDir->SetGuidance("Dose scorer commands");
0058 }
0059
0060
0061
0062 void Dose::SetNewValue(G4UIcommand* command, G4String newValue)
0063 {
0064 if (command == fpAddDoseCutOff.get()) {
0065 fDosesCutOff = G4UIcmdWithADoubleAndUnit::GetNewDoubleValue(newValue);
0066 }
0067 if (command == fpAddDoseToAbort.get()) {
0068 fDosesToAbort = G4UIcmdWithADoubleAndUnit::GetNewDoubleValue(newValue);
0069 }
0070 }
0071
0072
0073
0074 template<>
0075 void Scorer<Dose>::SetChemistryWorld(G4VChemistryWorld* pChemistryWorld)
0076 {
0077 fpChemistryWorld = pChemistryWorld;
0078 }
0079
0080
0081
0082 template<>
0083 G4VChemistryWorld* Scorer<Dose>::GetChemistryWorld() const
0084 {
0085 return fpChemistryWorld;
0086 }
0087
0088
0089
0090 template<>
0091 void Scorer<Dose>::clear()
0092 {
0093 fpScorer->fCumulatedDose = 0.;
0094 }
0095
0096
0097
0098 template<>
0099 void Scorer<Dose>::Initialize(G4HCofThisEvent* HCE)
0100 {
0101 clear();
0102 fpEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
0103 if (fHCID < 0) {
0104 fHCID = GetCollectionID(0);
0105 }
0106 HCE->AddHitsCollection(fHCID, (G4VHitsCollection*)fpEvtMap);
0107 }
0108
0109
0110
0111 template<>
0112 void Scorer<Dose>::EndOfEvent(G4HCofThisEvent*)
0113 {
0114 if (!G4RunManager::GetRunManager()->GetCurrentEvent()->IsAborted()) {
0115 fpEvtMap->add(0, fpScorer->fDosesCutOff);
0116 }
0117 fpScorer->fCumulatedDose = 0.;
0118 }
0119
0120
0121
0122 template<>
0123 G4bool Scorer<Dose>::ProcessHits(G4Step* aStep, G4TouchableHistory*)
0124 {
0125 auto currentEvent = G4EventManager::GetEventManager();
0126 const G4Track* track = aStep->GetTrack();
0127 auto boundingBox = fpChemistryWorld->GetChemistryBoundary();
0128 G4double V = boundingBox->Volume() / cm3;
0129 G4double edep = aStep->GetTotalEnergyDeposit();
0130 if (edep == 0.) {
0131 return false;
0132 }
0133 (fpScorer->fCumulatedDose) += edep;
0134 if (track->GetParentID() == 0 && aStep->IsFirstStepInVolume()) {
0135 G4double DoseInGray = ((fpScorer->fCumulatedDose) / eV) / (0.001 * V * 6.242e+18);
0136 if (DoseInGray > fpScorer->fDosesCutOff / gray) {
0137 G4cout << "_____________________________________________________________________________"
0138 << G4endl;
0139 auto name = currentEvent->GetConstCurrentEvent()
0140 ->GetPrimaryVertex()
0141 ->GetPrimary()
0142 ->GetParticleDefinition()
0143 ->GetParticleName();
0144 auto energy =
0145 currentEvent->GetConstCurrentEvent()->GetPrimaryVertex()->GetPrimary()->GetKineticEnergy();
0146 G4cout << "Stop this beam line (" << name << ", " << energy
0147 << " MeV) at actual dose: " << DoseInGray
0148 << " Gy. Cut-off dose: " << fpScorer->fDosesCutOff / gray << " Gy" << G4endl;
0149 G4cout << "The beam of " << 1000000 - currentEvent->GetStackManager()->GetNUrgentTrack() - 1
0150 << " tracks"
0151 << " in a volume of " << V * 1e+12
0152 << " um3. Total deposit energy: " << fpScorer->fCumulatedDose / eV << " eV. "
0153 << G4endl;
0154 if (DoseInGray > fpScorer->fDosesToAbort / gray) {
0155 G4cout << "Abort this beam line (" << name << ", " << energy
0156 << " MeV) at actual dose: " << DoseInGray << " Gy." << G4endl;
0157 G4RunManager::GetRunManager()->AbortEvent();
0158 }
0159 G4cout << "_____________________________________________________________________________"
0160 << G4endl;
0161 auto myTrack = ((G4Track*)track);
0162 myTrack->SetTrackStatus(fStopAndKill);
0163 auto secondaries = track->GetStep()->GetSecondaryInCurrentStep();
0164 if (!secondaries->empty()) {
0165 for (auto it : *(secondaries)) {
0166 if (it != nullptr) {
0167 ((G4Track*)it)->SetTrackStatus(fStopAndKill);
0168 }
0169 }
0170 }
0171 currentEvent->GetStackManager()->ClearUrgentStack();
0172 }
0173 }
0174 return true;
0175 }
0176
0177
0178
0179
0180 Gvalues::Gvalues()
0181 : G4UImessenger(),
0182 fTimeLimit(G4Scheduler::Instance()->GetEndTime()),
0183 fSpeciesdir(new G4UIdirectory("/scorer/Gvalues/")),
0184 fTimeBincmd(new G4UIcmdWithAnInteger("/scorer/Gvalues/nOfTimeBins", this)),
0185 fAddTimeToRecordcmd(new G4UIcmdWithADoubleAndUnit("/scorer/Gvalues/addTimeToRecord", this))
0186 {
0187 fSpeciesdir->SetGuidance("ScoreSpecies commands");
0188 G4MoleculeCounter::Instance()->ResetCounter();
0189 }
0190
0191
0192
0193 void Gvalues::SetNewValue(G4UIcommand* command, G4String newValue)
0194 {
0195 if (command == fAddTimeToRecordcmd.get()) {
0196 G4double cmdTime = fAddTimeToRecordcmd->GetNewDoubleValue(newValue);
0197 if (fTimeLimit >= cmdTime) {
0198 AddTimeToRecord(cmdTime);
0199 }
0200 else {
0201 AddTimeToRecord(fTimeLimit);
0202 }
0203 }
0204 if (command == fTimeBincmd.get()) {
0205 G4int cmdBins = fTimeBincmd->GetNewIntValue(newValue);
0206 G4double timeMin = 1 * ps;
0207 G4double timeMax = G4Scheduler::Instance()->GetEndTime() - timeMin;
0208 G4double timeLogMin = std::log10(timeMin);
0209 G4double timeLogMax = std::log10(timeMax);
0210 for (int i = 0; i <= cmdBins; i++) {
0211 AddTimeToRecord(std::pow(10, timeLogMin + i * (timeLogMax - timeLogMin) / cmdBins));
0212 }
0213 }
0214 }
0215
0216
0217
0218 void Gvalues::WriteWithAnalysisManager(G4VAnalysisManager* analysisManager, const std::string& out)
0219 {
0220 analysisManager->CreateNtuple(out, out);
0221 G4cout << "NtupleID : " << fRunID << " name : " << out << G4endl;
0222 analysisManager->CreateNtupleIColumn(fRunID, "speciesID");
0223 analysisManager->CreateNtupleIColumn(fRunID, "number");
0224 analysisManager->CreateNtupleIColumn(fRunID, "nEvent");
0225 analysisManager->CreateNtupleSColumn(fRunID, "speciesName");
0226 analysisManager->CreateNtupleDColumn(fRunID, "time");
0227 analysisManager->CreateNtupleDColumn(fRunID, "sumG");
0228 analysisManager->CreateNtupleDColumn(fRunID, "sumG2");
0229 analysisManager->FinishNtuple(fRunID);
0230
0231 for (const auto& it_map1 : fSpeciesInfoPerTime) {
0232 const InnerSpeciesMap& map2 = it_map1.second;
0233 for (auto it_map2 : map2) {
0234 double time = it_map1.first;
0235 auto species = it_map2.first;
0236 const G4String& name = species->GetName();
0237 int molID = it_map2.first->GetMoleculeID();
0238 auto number = it_map2.second.fNumber;
0239 double G = it_map2.second.fG;
0240 double G2 = it_map2.second.fG2;
0241
0242 analysisManager->FillNtupleIColumn(fRunID, 0, molID);
0243 analysisManager->FillNtupleIColumn(fRunID, 1, number);
0244 analysisManager->FillNtupleIColumn(fRunID, 2, fNEvent);
0245 analysisManager->FillNtupleSColumn(fRunID, 3, name);
0246 analysisManager->FillNtupleDColumn(fRunID, 4, time);
0247 analysisManager->FillNtupleDColumn(fRunID, 5, G);
0248 analysisManager->FillNtupleDColumn(fRunID, 6, G2);
0249 analysisManager->AddNtupleRow(fRunID);
0250 }
0251 }
0252 fRunID++;
0253 }
0254
0255
0256
0257 template<>
0258 void Scorer<Gvalues>::clear()
0259 {
0260 fpEvtMap->clear();
0261 fpScorer->fNEvent = 0;
0262 fpScorer->fEdep = 0;
0263 fpScorer->fSpeciesInfoPerTime.clear();
0264 }
0265
0266
0267
0268 template<>
0269 void Scorer<Gvalues>::Initialize(G4HCofThisEvent* HCE)
0270 {
0271 fpEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
0272 if (fHCID < 0) {
0273 fHCID = GetCollectionID(0);
0274 }
0275 HCE->AddHitsCollection(fHCID, (G4VHitsCollection*)fpEvtMap);
0276 }
0277
0278
0279
0280 template<>
0281 void Scorer<Gvalues>::SetChemistryWorld(G4VChemistryWorld* pChemistryWorld)
0282 {
0283 fpChemistryWorld = pChemistryWorld;
0284 }
0285
0286
0287
0288 template<>
0289 G4VChemistryWorld* Scorer<Gvalues>::GetChemistryWorld() const
0290 {
0291 return fpChemistryWorld;
0292 }
0293
0294
0295
0296 template<>
0297 G4bool Scorer<Gvalues>::ProcessHits(G4Step* aStep, G4TouchableHistory*)
0298 {
0299 G4double edep = aStep->GetTotalEnergyDeposit();
0300 if (edep == 0.) {
0301 return FALSE;
0302 }
0303 edep *= aStep->GetPreStepPoint()->GetWeight();
0304 G4int index = GetIndex(aStep);
0305 fpEvtMap->add(index, edep);
0306 (fpScorer->fEdep) += edep;
0307 return TRUE;
0308 }
0309
0310
0311
0312 template<>
0313 void Scorer<Gvalues>::SaveScavengerChange()
0314 {
0315 auto pScavengerMaterial =
0316 dynamic_cast<G4DNAScavengerMaterial*>(G4Scheduler::Instance()->GetScavengerMaterial());
0317 if (pScavengerMaterial == nullptr) {
0318 G4ExceptionDescription errMsg;
0319 errMsg << "pScavengerMaterial == nullptr";
0320 G4Exception("Scorer<Gvalues>::SaveScavengerChange()", "SaveScavengerChange",
0321 FatalErrorInArgument, errMsg);
0322 }
0323 auto scavengerList = pScavengerMaterial->GetScavengerList();
0324 auto V = fpChemistryWorld->GetChemistryBoundary()->Volume();
0325
0326 for (const auto& it : scavengerList) {
0327 if (it == G4MoleculeTable::Instance()->GetConfiguration("H2O")
0328 || G4MoleculeTable::Instance()->GetConfiguration("H3Op(B)") == it
0329 || G4MoleculeTable::Instance()->GetConfiguration("OHm(B)") == it)
0330 {
0331 continue;
0332 }
0333 for (auto time_mol : fpScorer->fTimeToRecord) {
0334 int64_t n_mol = pScavengerMaterial->GetNMoleculesAtTime(it, time_mol);
0335 if (n_mol < 0) {
0336 G4ExceptionDescription errMsg;
0337 errMsg << "SaveScavengerChange()::N molecules not valid < 0 : " << it->GetName()
0338 << " N : " << n_mol << G4endl;
0339 G4Exception("", "N<0", FatalException, errMsg);
0340 }
0341
0342 Gvalues::SpeciesInfo& molInfo = fpScorer->fSpeciesInfoPerTime[time_mol][it];
0343 molInfo.fNumber += n_mol;
0344 G4double gValue = n_mol / (Avogadro * V * 1.0e-6 );
0345 molInfo.fG += gValue;
0346 molInfo.fG2 += gValue * gValue;
0347 }
0348 }
0349 }
0350
0351
0352
0353 template<>
0354 void Scorer<Gvalues>::SaveMoleculeCounter()
0355 {
0356 if (fpEventScheduler == nullptr) {
0357 G4Exception("fpEventScheduler == nullptr", "Scorer<Gvalues>::SaveMoleculeCounter()",
0358 FatalException, "fpEventScheduler == nullptr");
0359 }
0360 else {
0361 auto counterMap = fpEventScheduler->GetCounterMap();
0362 if (counterMap.empty()) {
0363 if (!G4MoleculeCounter::Instance()->InUse()) {
0364 G4Exception("No counter", "Scorer<Gvalues>::SaveMoleculeCounter()", JustWarning,
0365 "G4MoleculeCounter::Instance() is not used");
0366 return;
0367 }
0368
0369 G4MoleculeCounter::RecordedMolecules species;
0370 species = G4MoleculeCounter::Instance()->GetRecordedMolecules();
0371 if (species.get() == nullptr) {
0372 return;
0373 }
0374 else if (species->empty()) {
0375 G4cout << "No molecule recorded, energy deposited" << G4endl;
0376 ++(fpScorer->fNEvent);
0377 fpScorer->fEdep = 0.;
0378 G4MoleculeCounter::Instance()->ResetCounter();
0379 return;
0380 }
0381 for (auto molecule : *species) {
0382 if (molecule == G4MoleculeTable::Instance()->GetConfiguration("O2")) {
0383 continue;
0384 }
0385 for (auto time_mol : fpScorer->fTimeToRecord) {
0386 int n_mol = G4MoleculeCounter::Instance()->GetNMoleculesAtTime(molecule, time_mol);
0387
0388 if (n_mol < 0) {
0389 G4ExceptionDescription errMsg;
0390 errMsg << "N molecules not valid < 0 " << G4endl;
0391 G4Exception("", "N<0", FatalException, errMsg);
0392 }
0393
0394 Gvalues::SpeciesInfo& molInfo = fpScorer->fSpeciesInfoPerTime[time_mol][molecule];
0395 molInfo.fNumber += n_mol;
0396 G4double gValue = (n_mol / (fpScorer->fEdep / eV)) * 100.;
0397 molInfo.fG += gValue;
0398 molInfo.fG2 += gValue * gValue;
0399 }
0400 }
0401 }
0402 else {
0403 for (const auto& map_mol : counterMap) {
0404 auto time_mol = map_mol.first;
0405 for (auto it_mol : map_mol.second) {
0406 auto molecule = it_mol.first;
0407 if (molecule == G4MoleculeTable::Instance()->GetConfiguration("O2")) {
0408 continue;
0409 }
0410 int n_mol = it_mol.second;
0411
0412 if (n_mol < 0) {
0413 G4ExceptionDescription errMsg;
0414 errMsg << "N molecules not valid < 0 "
0415 << " molecule : " << it_mol.first->GetName() << " N : " << n_mol << G4endl;
0416 G4Exception("", "N<0", FatalException, errMsg);
0417 }
0418
0419 Gvalues::SpeciesInfo& molInfo = fpScorer->fSpeciesInfoPerTime[time_mol][molecule];
0420 molInfo.fNumber += n_mol;
0421 G4double gValue = (n_mol / (fpScorer->fEdep / eV)) * 100.;
0422 molInfo.fG += gValue;
0423 molInfo.fG2 += gValue * gValue;
0424 }
0425 }
0426 }
0427 }
0428 }
0429
0430
0431
0432 template<>
0433 void Scorer<Gvalues>::EndOfEvent(G4HCofThisEvent*)
0434 {
0435 if (G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted()) {
0436 fpScorer->fEdep = 0.;
0437 G4MoleculeCounter::Instance()->ResetCounter();
0438 return;
0439 }
0440
0441 SaveScavengerChange();
0442 SaveMoleculeCounter();
0443
0444 ++(fpScorer->fNEvent);
0445 fpScorer->fEdep = 0.;
0446
0447 G4MoleculeCounter::Instance()->ResetCounter();
0448 G4MoleculeCounter::Instance()->Use(true);
0449 }
0450
0451
0452
0453 template<>
0454 void Scorer<Gvalues>::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
0455 {
0456 auto right = dynamic_cast<Scorer<Gvalues>*>(dynamic_cast<G4VPrimitiveScorer*>(workerScorer));
0457
0458 if (right == nullptr) {
0459 return;
0460 }
0461 if (right == this) {
0462 return;
0463 }
0464 auto it_map1 = right->fpScorer->fSpeciesInfoPerTime.begin();
0465 auto end_map1 = right->fpScorer->fSpeciesInfoPerTime.end();
0466
0467 for (; it_map1 != end_map1; ++it_map1) {
0468 Gvalues::InnerSpeciesMap& map2 = it_map1->second;
0469 auto it_map2 = map2.begin();
0470 auto end_map2 = map2.end();
0471
0472 for (; it_map2 != end_map2; ++it_map2) {
0473 Gvalues::SpeciesInfo& molInfo = fpScorer->fSpeciesInfoPerTime[it_map1->first][it_map2->first];
0474 molInfo.fNumber += it_map2->second.fNumber;
0475 molInfo.fG += it_map2->second.fG;
0476 molInfo.fG2 += it_map2->second.fG2;
0477 }
0478 }
0479 right->fpScorer->fSpeciesInfoPerTime.clear();
0480 fpScorer->fNEvent += right->fpScorer->fNEvent;
0481 right->fpScorer->fNEvent = 0;
0482 right->fpScorer->fEdep = 0.;
0483 }
0484
0485
0486
0487 template<>
0488 void Scorer<Gvalues>::OutputAndClear(const std::string& dose)
0489 {
0490 if (G4Threading::IsWorkerThread()) {
0491 return;
0492 }
0493 G4VAnalysisManager* analysisManager = G4AnalysisManager::Instance();
0494 if (analysisManager != nullptr) {
0495 this->fpScorer->WriteWithAnalysisManager(analysisManager, dose);
0496 }
0497 fpScorer->fNEvent = 0;
0498 fpScorer->fSpeciesInfoPerTime.clear();
0499 }
0500
0501